API reference

The public surface is small by design: one data model, one problem layer, one solve() front door.

Front door

Solve any QuGrid problem with any registered solver.

import qugrid as qg prob = qg.problems.Islanding(qg.cases.case9()) res = qg.solve(prob, solver="qaoa", seed=0) print(res.summary())

solver="auto" picks a sensible default: QAOA for small combinatorial problems, simulated annealing beyond statevector reach; HHL for small linear systems, exact algebra beyond.

repair="greedy" post-processes a combinatorial result with :func:qugrid.solvers.greedy_repair: the answer bitstring — and, for solvers that report an output distribution, every top state — is repaired toward feasibility, and P(optimum | repaired) lands in result.extras next to the raw P(optimum).

QUBO problems also accept the external stacks — "dimod-exact", "dwave-sa", "qiskit-qaoa" — which return the same :class:Result and need the matching extra installed.

Source code in src/qugrid/solvers/__init__.py

def solve(problem, solver: str = "auto", repair: str | None = None, **kwargs) -> Result:
    """Solve any QuGrid problem with any registered solver.

    >>> import qugrid as qg
    >>> prob = qg.problems.Islanding(qg.cases.case9())
    >>> res = qg.solve(prob, solver="qaoa", seed=0)
    >>> print(res.summary())

    ``solver="auto"`` picks a sensible default: QAOA for small combinatorial
    problems, simulated annealing beyond statevector reach; HHL for small
    linear systems, exact algebra beyond.

    ``repair="greedy"`` post-processes a combinatorial result with
    :func:`qugrid.solvers.greedy_repair`: the answer bitstring — and, for
    solvers that report an output distribution, every top state — is
    repaired toward feasibility, and ``P(optimum | repaired)`` lands in
    ``result.extras`` next to the raw ``P(optimum)``.

    QUBO problems also accept the external stacks — ``"dimod-exact"``,
    ``"dwave-sa"``, ``"qiskit-qaoa"`` — which return the same
    :class:`Result` and need the matching extra installed.
    """
    if isinstance(problem, CombinatorialProblem):
        kind = "qubo"
        if solver == "auto":
            solver = "qaoa" if problem.n <= 16 else "sa"
        if solver not in ("exact", "dimod-exact"):
            dr = problem.qubo.dynamic_range()
            if dr > 1e3:
                warnings.warn(
                    f"QUBO dynamic range is {dr:.2g} (>1e3): penalty terms dwarf the "
                    "cost differences, and sampling solvers may not resolve them. "
                    "Consider unbalanced penalties (QUBOBuilder.add_inequality) or "
                    "qugrid.methods.AugmentedLagrangianLoop.",
                    RuntimeWarning,
                    stacklevel=2,
                )
    elif isinstance(problem, LinearSystemProblem):
        kind = "linear"
        if solver == "auto":
            solver = "hhl" if problem.n <= 16 else "numpy"
    else:
        raise TypeError(
            f"solve() expects a CombinatorialProblem or LinearSystemProblem, "
            f"got {type(problem).__name__}"
        )
    if solver not in REGISTRY:
        options = ", ".join(name for name, (_, k) in REGISTRY.items() if k == kind)
        raise ValueError(f"unknown solver {solver!r}; options for this problem: {options}")
    fn, expected = REGISTRY[solver]
    if expected != kind:
        options = ", ".join(name for name, (_, k) in REGISTRY.items() if k == kind)
        raise ValueError(
            f"solver {solver!r} does not apply to {type(problem).__name__}; "
            f"options: {options}"
        )
    res = fn(problem, **kwargs)
    if repair is not None:
        if kind != "qubo":
            raise ValueError("repair applies to combinatorial problems only")
        if repair != "greedy":
            raise ValueError(f"unknown repair {repair!r}; the only option is 'greedy'")
        repair_result(res, problem)
    return res

Data model

A power network in MATPOWER case format with computed conveniences.

Source code in src/qugrid/network.py

@dataclass
class Network:
    """A power network in MATPOWER case format with computed conveniences."""

    baseMVA: float
    bus: np.ndarray
    gen: np.ndarray
    branch: np.ndarray
    gencost: np.ndarray | None = None
    name: str = "network"
    _ext2int: dict[int, int] = field(init=False, repr=False, default_factory=dict)

    def __post_init__(self) -> None:
        self.bus = np.atleast_2d(np.asarray(self.bus, dtype=float)).copy()
        self.gen = np.atleast_2d(np.asarray(self.gen, dtype=float)).copy()
        self.branch = np.atleast_2d(np.asarray(self.branch, dtype=float)).copy()
        if self.gencost is not None and np.size(self.gencost):
            self.gencost = np.atleast_2d(np.asarray(self.gencost, dtype=float)).copy()
        else:
            self.gencost = None
        self._ext2int = {int(b): i for i, b in enumerate(self.bus[:, idx.BUS_I])}
        if len(self._ext2int) != self.n_bus:
            raise ValueError("duplicate bus numbers in case data")

    # ------------------------------------------------------------------ sizes
    @property
    def n_bus(self) -> int:
        return self.bus.shape[0]

    @property
    def n_gen(self) -> int:
        return self.gen.shape[0]

    @property
    def n_branch(self) -> int:
        return self.branch.shape[0]

    # ------------------------------------------------------------- index maps
    def bus_index(self, external: int | np.ndarray) -> np.ndarray | int:
        """Map external bus number(s) to internal 0-based position(s)."""
        if np.isscalar(external):
            return self._ext2int[int(external)]
        return np.array([self._ext2int[int(b)] for b in np.asarray(external).ravel()])

    @property
    def f_bus(self) -> np.ndarray:
        """Internal index of each branch's from-bus."""
        return self.bus_index(self.branch[:, idx.F_BUS])

    @property
    def t_bus(self) -> np.ndarray:
        """Internal index of each branch's to-bus."""
        return self.bus_index(self.branch[:, idx.T_BUS])

    @property
    def gen_bus(self) -> np.ndarray:
        """Internal bus index of each generator."""
        return self.bus_index(self.gen[:, idx.GEN_BUS])

    @property
    def branch_on(self) -> np.ndarray:
        return self.branch[:, idx.BR_STATUS] > 0

    @property
    def gen_on(self) -> np.ndarray:
        return self.gen[:, idx.GEN_STATUS] > 0

    @property
    def ref(self) -> int:
        """Internal index of the reference (slack) bus."""
        refs = np.flatnonzero(self.bus[:, idx.BUS_TYPE] == idx.REF)
        if len(refs) == 0:
            raise ValueError("case has no reference bus")
        return int(refs[0])

    @property
    def pv(self) -> np.ndarray:
        return np.flatnonzero(self.bus[:, idx.BUS_TYPE] == idx.PV)

    @property
    def pq(self) -> np.ndarray:
        return np.flatnonzero(self.bus[:, idx.BUS_TYPE] == idx.PQ)

    # ------------------------------------------------------------ injections
    @property
    def load_p(self) -> np.ndarray:
        """Active load per bus [MW]."""
        return self.bus[:, idx.PD].copy()

    @property
    def load_q(self) -> np.ndarray:
        """Reactive load per bus [MVAr]."""
        return self.bus[:, idx.QD].copy()

    def gen_p_per_bus(self) -> np.ndarray:
        """Scheduled in-service generation per bus [MW]."""
        p = np.zeros(self.n_bus)
        on = self.gen_on
        np.add.at(p, self.gen_bus[on], self.gen[on, idx.PG])
        return p

    def sbus(self) -> np.ndarray:
        """Complex net power injection per bus [pu]."""
        sload = (self.bus[:, idx.PD] + 1j * self.bus[:, idx.QD]) / self.baseMVA
        sgen = np.zeros(self.n_bus, dtype=complex)
        on = self.gen_on
        np.add.at(
            sgen,
            self.gen_bus[on],
            (self.gen[on, idx.PG] + 1j * self.gen[on, idx.QG]) / self.baseMVA,
        )
        return sgen - sload

    # ---------------------------------------------------------------- Y / Bdc
    def ybus(self) -> np.ndarray:
        """Dense complex bus admittance matrix [pu].

        Follows MATPOWER's ``makeYbus`` (branch pi model with off-nominal tap
        ratio and phase shift). Dense is a deliberate choice: QuGrid targets
        research-scale cases (up to a few hundred buses), where dense algebra
        is simpler and fast enough.
        """
        nb, nl = self.n_bus, self.n_branch
        stat = self.branch_on.astype(float)
        ys = stat / (self.branch[:, idx.BR_R] + 1j * self.branch[:, idx.BR_X])
        bc = stat * self.branch[:, idx.BR_B]
        tap_mag = np.where(self.branch[:, idx.TAP] == 0.0, 1.0, self.branch[:, idx.TAP])
        tap = tap_mag * np.exp(1j * np.deg2rad(self.branch[:, idx.SHIFT]))
        ytt = ys + 1j * bc / 2
        yff = ytt / (tap * np.conj(tap))
        yft = -ys / np.conj(tap)
        ytf = -ys / tap
        ysh = (self.bus[:, idx.GS] + 1j * self.bus[:, idx.BS]) / self.baseMVA
        f, t = self.f_bus, self.t_bus
        y = np.zeros((nb, nb), dtype=complex)
        np.add.at(y, (f, f), yff)
        np.add.at(y, (t, t), ytt)
        np.add.at(y, (f, t), yft)
        np.add.at(y, (t, f), ytf)
        y[np.arange(nb), np.arange(nb)] += ysh
        del nl
        return y

    def bdc(self) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
        """DC power flow matrices, following MATPOWER's ``makeBdc``.

        Returns ``(Bbus, Bf, Pbusinj, Pfinj)`` such that, in per unit,
        ``P_bus = Bbus @ theta + Pbusinj`` and ``P_branch = Bf @ theta + Pfinj``.
        """
        nb = self.n_bus
        stat = self.branch_on.astype(float)
        b = stat / self.branch[:, idx.BR_X]
        tap_mag = np.where(self.branch[:, idx.TAP] == 0.0, 1.0, self.branch[:, idx.TAP])
        b = b / tap_mag
        f, t = self.f_bus, self.t_bus
        nl = self.n_branch
        bf = np.zeros((nl, nb))
        bf[np.arange(nl), f] += b
        bf[np.arange(nl), t] -= b
        bbus = np.zeros((nb, nb))
        np.add.at(bbus, (f, f), b)
        np.add.at(bbus, (t, t), b)
        np.add.at(bbus, (f, t), -b)
        np.add.at(bbus, (t, f), -b)
        pfinj = b * (-np.deg2rad(self.branch[:, idx.SHIFT]))
        pbusinj = np.zeros(nb)
        np.add.at(pbusinj, f, pfinj)
        np.add.at(pbusinj, t, -pfinj)
        return bbus, bf, pbusinj, pfinj

    # ------------------------------------------------------------------ graph
    def edges(self, in_service_only: bool = True) -> list[tuple[int, int, float]]:
        """Branch list as ``(from, to, |b|)`` with internal indices.

        The weight is the branch susceptance magnitude ``1/x`` — the standard
        measure of electrical coupling used in partitioning studies.
        """
        out = []
        for line in range(self.n_branch):
            if in_service_only and not self.branch_on[line]:
                continue
            w = abs(1.0 / self.branch[line, idx.BR_X])
            out.append((int(self.f_bus[line]), int(self.t_bus[line]), float(w)))
        return out

    def adjacency(self) -> np.ndarray:
        """Symmetric 0/1 adjacency matrix of the in-service network."""
        a = np.zeros((self.n_bus, self.n_bus))
        for f, t, _ in self.edges():
            a[f, t] = a[t, f] = 1.0
        return a

    def is_connected(self) -> bool:
        a = self.adjacency()
        seen = np.zeros(self.n_bus, dtype=bool)
        stack = [0]
        seen[0] = True
        while stack:
            u = stack.pop()
            for v in np.flatnonzero(a[u]):
                if not seen[v]:
                    seen[v] = True
                    stack.append(int(v))
        return bool(seen.all())

    # ------------------------------------------------------------- transforms
    def copy(self) -> Network:
        return Network(
            baseMVA=self.baseMVA,
            bus=self.bus,
            gen=self.gen,
            branch=self.branch,
            gencost=self.gencost,
            name=self.name,
        )

    def scale_loads(self, factor: float | np.ndarray) -> Network:
        """Return a copy with active and reactive loads multiplied by ``factor``.

        ``factor`` is a scalar or a per-bus array of length ``n_bus``.
        """
        out = self.copy()
        out.bus[:, idx.PD] *= factor
        out.bus[:, idx.QD] *= factor
        return out

    def drop_branch(self, line: int) -> Network:
        """Return a copy with one branch switched out (N-1 outage)."""
        out = self.copy()
        out.branch[line, idx.BR_STATUS] = 0
        return out

    # ------------------------------------------------------------ constructors
    @classmethod
    def from_ppc(cls, ppc: dict, name: str = "network") -> Network:
        """Build from a MATPOWER/PYPOWER case dict (``baseMVA, bus, gen, branch``)."""
        return cls(
            baseMVA=float(ppc["baseMVA"]),
            bus=ppc["bus"],
            gen=ppc["gen"],
            branch=ppc["branch"],
            gencost=ppc.get("gencost"),
            name=name,
        )

    @classmethod
    def from_matpower(cls, path: str, name: str | None = None) -> Network:
        """Read a MATPOWER ``.m`` case file."""
        from pathlib import Path

        from qugrid.io.matpower import read_matpower

        ppc = read_matpower(path)
        return cls.from_ppc(ppc, name=name or Path(path).stem)

    @classmethod
    def from_pandapower(cls, net, name: str | None = None) -> Network:
        """Convert a pandapower network (requires ``pandapower`` installed)."""
        try:
            from pandapower.converter import to_ppc
        except ImportError as err:  # pragma: no cover - exercised without extra
            raise ImportError(
                "pandapower is required for Network.from_pandapower; "
                "install with `pip install qugrid[pandapower]`"
            ) from err
        ppc = to_ppc(net, init="flat")
        return cls.from_ppc(ppc, name=name or "pandapower-net")

    def to_ppc(self) -> dict:
        """Export as a PYPOWER-compatible case dict."""
        ppc = {
            "version": "2",
            "baseMVA": self.baseMVA,
            "bus": self.bus.copy(),
            "gen": self.gen.copy(),
            "branch": self.branch.copy(),
        }
        if self.gencost is not None:
            ppc["gencost"] = self.gencost.copy()
        return ppc

    # ------------------------------------------------------------------ report
    def __repr__(self) -> str:
        return (
            f"Network({self.name!r}: {self.n_bus} buses, {self.n_branch} branches, "
            f"{self.n_gen} generators, load {self.load_p.sum():.1f} MW)"
        )

f_bus property

Internal index of each branch's from-bus.

t_bus property

Internal index of each branch's to-bus.

gen_bus property

Internal bus index of each generator.

ref property

Internal index of the reference (slack) bus.

load_p property

Active load per bus [MW].

load_q property

Reactive load per bus [MVAr].

bus_index(external)

Map external bus number(s) to internal 0-based position(s).

Source code in src/qugrid/network.py

def bus_index(self, external: int | np.ndarray) -> np.ndarray | int:
    """Map external bus number(s) to internal 0-based position(s)."""
    if np.isscalar(external):
        return self._ext2int[int(external)]
    return np.array([self._ext2int[int(b)] for b in np.asarray(external).ravel()])

gen_p_per_bus()

Scheduled in-service generation per bus [MW].

Source code in src/qugrid/network.py

def gen_p_per_bus(self) -> np.ndarray:
    """Scheduled in-service generation per bus [MW]."""
    p = np.zeros(self.n_bus)
    on = self.gen_on
    np.add.at(p, self.gen_bus[on], self.gen[on, idx.PG])
    return p

sbus()

Complex net power injection per bus [pu].

Source code in src/qugrid/network.py

def sbus(self) -> np.ndarray:
    """Complex net power injection per bus [pu]."""
    sload = (self.bus[:, idx.PD] + 1j * self.bus[:, idx.QD]) / self.baseMVA
    sgen = np.zeros(self.n_bus, dtype=complex)
    on = self.gen_on
    np.add.at(
        sgen,
        self.gen_bus[on],
        (self.gen[on, idx.PG] + 1j * self.gen[on, idx.QG]) / self.baseMVA,
    )
    return sgen - sload

ybus()

Dense complex bus admittance matrix [pu].

Follows MATPOWER's makeYbus (branch pi model with off-nominal tap ratio and phase shift). Dense is a deliberate choice: QuGrid targets research-scale cases (up to a few hundred buses), where dense algebra is simpler and fast enough.

Source code in src/qugrid/network.py

def ybus(self) -> np.ndarray:
    """Dense complex bus admittance matrix [pu].

    Follows MATPOWER's ``makeYbus`` (branch pi model with off-nominal tap
    ratio and phase shift). Dense is a deliberate choice: QuGrid targets
    research-scale cases (up to a few hundred buses), where dense algebra
    is simpler and fast enough.
    """
    nb, nl = self.n_bus, self.n_branch
    stat = self.branch_on.astype(float)
    ys = stat / (self.branch[:, idx.BR_R] + 1j * self.branch[:, idx.BR_X])
    bc = stat * self.branch[:, idx.BR_B]
    tap_mag = np.where(self.branch[:, idx.TAP] == 0.0, 1.0, self.branch[:, idx.TAP])
    tap = tap_mag * np.exp(1j * np.deg2rad(self.branch[:, idx.SHIFT]))
    ytt = ys + 1j * bc / 2
    yff = ytt / (tap * np.conj(tap))
    yft = -ys / np.conj(tap)
    ytf = -ys / tap
    ysh = (self.bus[:, idx.GS] + 1j * self.bus[:, idx.BS]) / self.baseMVA
    f, t = self.f_bus, self.t_bus
    y = np.zeros((nb, nb), dtype=complex)
    np.add.at(y, (f, f), yff)
    np.add.at(y, (t, t), ytt)
    np.add.at(y, (f, t), yft)
    np.add.at(y, (t, f), ytf)
    y[np.arange(nb), np.arange(nb)] += ysh
    del nl
    return y

bdc()

DC power flow matrices, following MATPOWER's makeBdc.

Returns (Bbus, Bf, Pbusinj, Pfinj) such that, in per unit, P_bus = Bbus @ theta + Pbusinj and P_branch = Bf @ theta + Pfinj.

Source code in src/qugrid/network.py

def bdc(self) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
    """DC power flow matrices, following MATPOWER's ``makeBdc``.

    Returns ``(Bbus, Bf, Pbusinj, Pfinj)`` such that, in per unit,
    ``P_bus = Bbus @ theta + Pbusinj`` and ``P_branch = Bf @ theta + Pfinj``.
    """
    nb = self.n_bus
    stat = self.branch_on.astype(float)
    b = stat / self.branch[:, idx.BR_X]
    tap_mag = np.where(self.branch[:, idx.TAP] == 0.0, 1.0, self.branch[:, idx.TAP])
    b = b / tap_mag
    f, t = self.f_bus, self.t_bus
    nl = self.n_branch
    bf = np.zeros((nl, nb))
    bf[np.arange(nl), f] += b
    bf[np.arange(nl), t] -= b
    bbus = np.zeros((nb, nb))
    np.add.at(bbus, (f, f), b)
    np.add.at(bbus, (t, t), b)
    np.add.at(bbus, (f, t), -b)
    np.add.at(bbus, (t, f), -b)
    pfinj = b * (-np.deg2rad(self.branch[:, idx.SHIFT]))
    pbusinj = np.zeros(nb)
    np.add.at(pbusinj, f, pfinj)
    np.add.at(pbusinj, t, -pfinj)
    return bbus, bf, pbusinj, pfinj

edges(in_service_only=True)

Branch list as (from, to, |b|) with internal indices.

The weight is the branch susceptance magnitude 1/x — the standard measure of electrical coupling used in partitioning studies.

Source code in src/qugrid/network.py

def edges(self, in_service_only: bool = True) -> list[tuple[int, int, float]]:
    """Branch list as ``(from, to, |b|)`` with internal indices.

    The weight is the branch susceptance magnitude ``1/x`` — the standard
    measure of electrical coupling used in partitioning studies.
    """
    out = []
    for line in range(self.n_branch):
        if in_service_only and not self.branch_on[line]:
            continue
        w = abs(1.0 / self.branch[line, idx.BR_X])
        out.append((int(self.f_bus[line]), int(self.t_bus[line]), float(w)))
    return out

adjacency()

Symmetric 0/1 adjacency matrix of the in-service network.

Source code in src/qugrid/network.py

def adjacency(self) -> np.ndarray:
    """Symmetric 0/1 adjacency matrix of the in-service network."""
    a = np.zeros((self.n_bus, self.n_bus))
    for f, t, _ in self.edges():
        a[f, t] = a[t, f] = 1.0
    return a

scale_loads(factor)

Return a copy with active and reactive loads multiplied by factor.

factor is a scalar or a per-bus array of length n_bus.

Source code in src/qugrid/network.py

def scale_loads(self, factor: float | np.ndarray) -> Network:
    """Return a copy with active and reactive loads multiplied by ``factor``.

    ``factor`` is a scalar or a per-bus array of length ``n_bus``.
    """
    out = self.copy()
    out.bus[:, idx.PD] *= factor
    out.bus[:, idx.QD] *= factor
    return out

drop_branch(line)

Return a copy with one branch switched out (N-1 outage).

Source code in src/qugrid/network.py

def drop_branch(self, line: int) -> Network:
    """Return a copy with one branch switched out (N-1 outage)."""
    out = self.copy()
    out.branch[line, idx.BR_STATUS] = 0
    return out

from_ppc(ppc, name='network') classmethod

Build from a MATPOWER/PYPOWER case dict (baseMVA, bus, gen, branch).

Source code in src/qugrid/network.py

@classmethod
def from_ppc(cls, ppc: dict, name: str = "network") -> Network:
    """Build from a MATPOWER/PYPOWER case dict (``baseMVA, bus, gen, branch``)."""
    return cls(
        baseMVA=float(ppc["baseMVA"]),
        bus=ppc["bus"],
        gen=ppc["gen"],
        branch=ppc["branch"],
        gencost=ppc.get("gencost"),
        name=name,
    )

from_matpower(path, name=None) classmethod

Read a MATPOWER .m case file.

Source code in src/qugrid/network.py

@classmethod
def from_matpower(cls, path: str, name: str | None = None) -> Network:
    """Read a MATPOWER ``.m`` case file."""
    from pathlib import Path

    from qugrid.io.matpower import read_matpower

    ppc = read_matpower(path)
    return cls.from_ppc(ppc, name=name or Path(path).stem)

from_pandapower(net, name=None) classmethod

Convert a pandapower network (requires pandapower installed).

Source code in src/qugrid/network.py

@classmethod
def from_pandapower(cls, net, name: str | None = None) -> Network:
    """Convert a pandapower network (requires ``pandapower`` installed)."""
    try:
        from pandapower.converter import to_ppc
    except ImportError as err:  # pragma: no cover - exercised without extra
        raise ImportError(
            "pandapower is required for Network.from_pandapower; "
            "install with `pip install qugrid[pandapower]`"
        ) from err
    ppc = to_ppc(net, init="flat")
    return cls.from_ppc(ppc, name=name or "pandapower-net")

to_ppc()

Export as a PYPOWER-compatible case dict.

Source code in src/qugrid/network.py

def to_ppc(self) -> dict:
    """Export as a PYPOWER-compatible case dict."""
    ppc = {
        "version": "2",
        "baseMVA": self.baseMVA,
        "bus": self.bus.copy(),
        "gen": self.gen.copy(),
        "branch": self.branch.copy(),
    }
    if self.gencost is not None:
        ppc["gencost"] = self.gencost.copy()
    return ppc

Problem formulations

Bases: CombinatorialProblem

Multi-period unit commitment with binary-encoded dispatch.

Parameters

gens: Generator parameters (cost coefficients in $/h, limits in MW). demand: Demand per period [MW]. power_bits: Bits per unit and period for the dispatch expansion. 2 gives four output levels between pmin and pmax. weight_balance, weight_gate: Penalty weights. None picks values that dominate the worst-case cost range, which is sufficient for exactness on feasible instances. initial_on: Commitment state before the first period (for startup costs).

Source code in src/qugrid/problems/unit_commitment.py

class UnitCommitment(CombinatorialProblem):
    """Multi-period unit commitment with binary-encoded dispatch.

    Parameters
    ----------
    gens:
        Generator parameters (cost coefficients in $/h, limits in MW).
    demand:
        Demand per period [MW].
    power_bits:
        Bits per unit and period for the dispatch expansion. ``2`` gives four
        output levels between ``pmin`` and ``pmax``.
    weight_balance, weight_gate:
        Penalty weights. ``None`` picks values that dominate the worst-case
        cost range, which is sufficient for exactness on feasible instances.
    initial_on:
        Commitment state before the first period (for startup costs).
    """

    def __init__(
        self,
        gens: list[GenParams],
        demand: list[float] | np.ndarray,
        power_bits: int = 2,
        weight_balance: float | None = None,
        weight_gate: float | None = None,
        initial_on: np.ndarray | None = None,
    ) -> None:
        self.gens = list(gens)
        self.demand = np.asarray(demand, dtype=float)
        self.power_bits = int(power_bits)
        if self.power_bits < 1:
            raise ValueError("power_bits must be >= 1")
        self.n_gen = len(self.gens)
        self.n_t = len(self.demand)
        self.initial_on = (
            np.zeros(self.n_gen) if initial_on is None else np.asarray(initial_on, float)
        )
        self.delta = np.array(
            [(g.pmax - g.pmin) / (2**self.power_bits - 1) for g in self.gens]
        )
        cost_ceiling = float(
            sum(g.cost(g.pmax) + g.startup for g in self.gens) * self.n_t
        )
        min_quantum = float(self.delta.min())
        self.weight_balance = (
            weight_balance
            if weight_balance is not None
            else 2.0 * cost_ceiling / max(min_quantum**2, 1e-9)
        )
        self.weight_gate = weight_gate if weight_gate is not None else 2.0 * cost_ceiling
        self.balance_tol = 0.75 * min_quantum

    # ------------------------------------------------------------- variables
    def _u(self, g: int, t: int) -> str:
        return f"u[{g},{t}]"

    def _b(self, g: int, t: int, k: int) -> str:
        return f"b[{g},{t},{k}]"

    @cached_property
    def qubo(self) -> QUBO:
        bld = QUBOBuilder()
        u = {(g, t): bld.var(self._u(g, t)) for g in range(self.n_gen) for t in range(self.n_t)}
        b = {
            (g, t, k): bld.var(self._b(g, t, k))
            for g in range(self.n_gen)
            for t in range(self.n_t)
            for k in range(self.power_bits)
        }

        for t in range(self.n_t):
            balance_terms: list[tuple[int, float]] = []
            for g, gen in enumerate(self.gens):
                # P_{g,t} as linear terms over (u, bits)
                p_terms = [(u[g, t], gen.pmin)] + [
                    (b[g, t, k], self.delta[g] * 2**k) for k in range(self.power_bits)
                ]
                balance_terms.extend(p_terms)
                # cost: c2 * P^2 + c1 * P + c0 * u
                bld.add_squared_penalty(p_terms, 0.0, weight=gen.c2)
                for i, coef in p_terms:
                    bld.add_linear(i, gen.c1 * coef)
                bld.add_linear(u[g, t], gen.c0)
                # startup: su * u_t * (1 - u_{t-1})
                if gen.startup:
                    if t == 0:
                        bld.add_linear(u[g, t], gen.startup * (1.0 - self.initial_on[g]))
                    else:
                        bld.add_linear(u[g, t], gen.startup)
                        bld.add_quadratic(u[g, t], u[g, t - 1], -gen.startup)
                # gating: bits force u on
                for k in range(self.power_bits):
                    bld.add_linear(b[g, t, k], self.weight_gate)
                    bld.add_quadratic(b[g, t, k], u[g, t], -self.weight_gate)
            bld.add_squared_penalty(balance_terms, -float(self.demand[t]), self.weight_balance)
        return bld.build()

    # --------------------------------------------------------------- decoding
    def _split(self, x: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
        names = self.qubo.names
        commit = np.zeros((self.n_gen, self.n_t))
        power = np.zeros((self.n_gen, self.n_t))
        val = {name: int(x[i]) for i, name in enumerate(names)}
        for g in range(self.n_gen):
            for t in range(self.n_t):
                commit[g, t] = val[self._u(g, t)]
                level = sum(2**k * val[self._b(g, t, k)] for k in range(self.power_bits))
                power[g, t] = self.gens[g].pmin * commit[g, t] + self.delta[g] * level
        return commit, power

    def decode(self, x: np.ndarray) -> dict:
        commit, power = self._split(x)
        cost = 0.0
        for t in range(self.n_t):
            for g, gen in enumerate(self.gens):
                if commit[g, t]:
                    cost += float(gen.cost(power[g, t]))
                prev = self.initial_on[g] if t == 0 else commit[g, t - 1]
                cost += gen.startup * max(commit[g, t] - prev, 0.0)
        balance = power.sum(axis=0) - self.demand
        return {
            "commit": commit,
            "power": power,
            "cost": cost,
            "balance_error_mw": float(np.abs(balance).max()),
        }

    def is_feasible(self, x: np.ndarray) -> bool:
        commit, power = self._split(x)
        gate_ok = bool(np.all(power[commit == 0] == 0))
        balance = power.sum(axis=0) - self.demand
        return gate_ok and bool(np.abs(balance).max() <= self.balance_tol)

    def constraint_residual(self, x: np.ndarray) -> float:
        """MW of violation: imbalance beyond tolerance plus gated-off power."""
        commit, power = self._split(x)
        balance = np.abs(power.sum(axis=0) - self.demand)
        gate = float(power[commit == 0].sum())
        return float(np.maximum(balance - self.balance_tol, 0.0).sum() + gate)

    # -------------------------------------------------------------- baselines
    def continuous_reference(self) -> dict:
        """True UC optimum with continuous dispatch (enumeration + exact ED)."""
        sol = solve_uc_enumerate(self.gens, self.demand, initial_on=self.initial_on)
        return {"commit": sol.commit, "power": sol.power, "cost": sol.cost}

constraint_residual(x)

MW of violation: imbalance beyond tolerance plus gated-off power.

Source code in src/qugrid/problems/unit_commitment.py

def constraint_residual(self, x: np.ndarray) -> float:
    """MW of violation: imbalance beyond tolerance plus gated-off power."""
    commit, power = self._split(x)
    balance = np.abs(power.sum(axis=0) - self.demand)
    gate = float(power[commit == 0].sum())
    return float(np.maximum(balance - self.balance_tol, 0.0).sum() + gate)

continuous_reference()

True UC optimum with continuous dispatch (enumeration + exact ED).

Source code in src/qugrid/problems/unit_commitment.py

def continuous_reference(self) -> dict:
    """True UC optimum with continuous dispatch (enumeration + exact ED)."""
    sol = solve_uc_enumerate(self.gens, self.demand, initial_on=self.initial_on)
    return {"commit": sol.commit, "power": sol.power, "cost": sol.cost}

Bases: CombinatorialProblem

Minimize generation cost subject to a demand-balance penalty.

Output of unit g: P_g = pmin_g + delta_g * sum_k 2^k b_{g,k} with delta_g = (pmax - pmin) / (2^K - 1). All units are committed.

Source code in src/qugrid/problems/economic_dispatch.py

class EconomicDispatchQUBO(CombinatorialProblem):
    """Minimize generation cost subject to a demand-balance penalty.

    Output of unit ``g``: ``P_g = pmin_g + delta_g * sum_k 2^k b_{g,k}`` with
    ``delta_g = (pmax - pmin) / (2^K - 1)``. All units are committed.
    """

    def __init__(
        self,
        gens: list[GenParams],
        demand: float,
        power_bits: int = 3,
        weight_balance: float | None = None,
    ) -> None:
        self.gens = list(gens)
        self.demand = float(demand)
        self.power_bits = int(power_bits)
        self.delta = np.array(
            [(g.pmax - g.pmin) / (2**self.power_bits - 1) for g in self.gens]
        )
        cost_ceiling = float(sum(g.cost(g.pmax) for g in self.gens))
        min_quantum = float(self.delta.min())
        self.weight_balance = (
            weight_balance
            if weight_balance is not None
            else 2.0 * cost_ceiling / max(min_quantum**2, 1e-9)
        )
        self.balance_tol = 0.75 * min_quantum

    @cached_property
    def qubo(self) -> QUBO:
        bld = QUBOBuilder()
        balance_terms: list[tuple[int, float]] = []
        offset = 0.0
        for g, gen in enumerate(self.gens):
            bits = [bld.var(f"b[{g},{k}]") for k in range(self.power_bits)]
            p_terms = [(bits[k], self.delta[g] * 2**k) for k in range(self.power_bits)]
            # cost of (pmin + sum) : expand c2*(pmin + s)^2 + c1*(pmin + s) + c0
            bld.add_squared_penalty(p_terms, gen.pmin, weight=gen.c2)
            for i, coef in p_terms:
                bld.add_linear(i, gen.c1 * coef)
            offset += gen.c1 * gen.pmin + gen.c0  # constant cost terms (all units on)
            balance_terms.extend(p_terms)
        total_pmin = float(sum(g.pmin for g in self.gens))
        bld.add_squared_penalty(
            balance_terms, total_pmin - self.demand, weight=self.weight_balance
        )
        bld.add_constant(offset)
        return bld.build()

    def power(self, x: np.ndarray) -> np.ndarray:
        x = np.asarray(x, dtype=int)
        p = np.zeros(len(self.gens))
        for g, gen in enumerate(self.gens):
            level = sum(
                2**k * x[g * self.power_bits + k] for k in range(self.power_bits)
            )
            p[g] = gen.pmin + self.delta[g] * level
        return p

    def decode(self, x: np.ndarray) -> dict:
        p = self.power(x)
        cost = float(sum(g.cost(p[i]) for i, g in enumerate(self.gens)))
        return {
            "power": p,
            "cost": cost,
            "balance_error_mw": float(abs(p.sum() - self.demand)),
        }

    def is_feasible(self, x: np.ndarray) -> bool:
        return bool(abs(self.power(x).sum() - self.demand) <= self.balance_tol)

    def continuous_reference(self) -> dict:
        p, cost = economic_dispatch(self.gens, self.demand)
        return {"power": p, "cost": cost}

Bases: CombinatorialProblem

Two-way controlled islanding of a network.

Variables: x_i = 0/1 assigns bus i to island A or B.

Parameters

net: The network. One binary variable per bus. alpha: Weight of the squared power-imbalance term (per-unit power). None scales it so a one-per-unit imbalance costs as much as cutting the strongest line. beta: Weight of the size-balance term. None uses a small default that breaks the all-one-island degeneracy without dominating.

.. note:: The no-split assignment (every bus in one island) is always a QUBO state; on very small or tightly meshed networks it can be the ground state under the default weights, because cutting any line costs more than the small size-balance reward. That optimum is infeasible for islanding — is_feasible is False and decode()["sizes"] shows an empty island. On the bundled IEEE cases (9 buses and larger) the defaults give a genuine split; if your network comes back trivial, raise beta (1.0 x mean edge weight forces a split on every bundled case) or fix seed buses by editing the QUBO.

Source code in src/qugrid/problems/islanding.py

class Islanding(CombinatorialProblem):
    """Two-way controlled islanding of a network.

    Variables: ``x_i = 0/1`` assigns bus ``i`` to island A or B.

    Parameters
    ----------
    net:
        The network. One binary variable per bus.
    alpha:
        Weight of the squared power-imbalance term (per-unit power). ``None``
        scales it so a one-per-unit imbalance costs as much as cutting the
        strongest line.
    beta:
        Weight of the size-balance term. ``None`` uses a small default that
        breaks the all-one-island degeneracy without dominating.

    .. note::
        The no-split assignment (every bus in one island) is always a QUBO
        state; on very small or tightly meshed networks it can be the ground
        state under the default weights, because cutting any line costs more
        than the small size-balance reward. That optimum is *infeasible* for
        islanding — ``is_feasible`` is False and ``decode()["sizes"]`` shows
        an empty island. On the bundled IEEE cases (9 buses and larger) the
        defaults give a genuine split; if your network comes back trivial,
        raise ``beta`` (1.0 x mean edge weight forces a split on every
        bundled case) or fix seed buses by editing the QUBO.
    """

    def __init__(self, net: Network, alpha: float | None = None, beta: float | None = None):
        self.net = net
        self.edge_list = net.edges()
        w = np.array([w for _, _, w in self.edge_list])
        self.injection_pu = (net.gen_p_per_bus() - net.load_p) / net.baseMVA
        self.alpha = alpha if alpha is not None else float(2.0 * w.max())
        self.beta = beta if beta is not None else float(0.05 * w.mean())

    @cached_property
    def qubo(self) -> QUBO:
        n = self.net.n_bus
        bld = QUBOBuilder()
        x = [bld.var(f"bus{int(self.net.bus[i, 0])}") for i in range(n)]
        # cut term: w * (x_f + x_t - 2 x_f x_t)
        for f, t, w in self.edge_list:
            bld.add_linear(x[f], w)
            bld.add_linear(x[t], w)
            bld.add_quadratic(x[f], x[t], -2.0 * w)
        # power imbalance: alpha * (sum_i p_i * s_i)^2 with s = 2x - 1
        terms = [(x[i], 2.0 * self.injection_pu[i]) for i in range(n)]
        const = -float(self.injection_pu.sum())
        bld.add_squared_penalty(terms, const, weight=self.alpha)
        # size balance: beta * (sum_i s_i)^2
        bld.add_squared_penalty([(x[i], 2.0) for i in range(n)], -float(n), weight=self.beta)
        return bld.build()

    # ---------------------------------------------------------------- decode
    def decode(self, x: np.ndarray) -> dict:
        x = np.asarray(x, dtype=int)
        cut = [(f, t) for f, t, _ in self.edge_list if x[f] != x[t]]
        cut_weight = float(sum(w for f, t, w in self.edge_list if x[f] != x[t]))
        p_a = float(self.injection_pu[x == 0].sum() * self.net.baseMVA)
        p_b = float(self.injection_pu[x == 1].sum() * self.net.baseMVA)
        return {
            "islands": x,
            "sizes": (int((x == 0).sum()), int((x == 1).sum())),
            "cut_lines": cut,
            "n_cut": len(cut),
            "cut_weight": cut_weight,
            "island_power_mw": (p_a, p_b),
            "islands_connected": self._connected(x),
        }

    def _connected(self, x: np.ndarray) -> tuple[bool, bool]:
        out = []
        for side in (0, 1):
            nodes = set(np.flatnonzero(x == side).tolist())
            if not nodes:
                out.append(False)
                continue
            adj: dict[int, list[int]] = {v: [] for v in nodes}
            for f, t, _ in self.edge_list:
                if f in nodes and t in nodes:
                    adj[f].append(t)
                    adj[t].append(f)
            start = next(iter(nodes))
            seen = {start}
            stack = [start]
            while stack:
                u = stack.pop()
                for v in adj[u]:
                    if v not in seen:
                        seen.add(v)
                        stack.append(v)
            out.append(seen == nodes)
        return (out[0], out[1])

    def is_feasible(self, x: np.ndarray) -> bool:
        x = np.asarray(x, dtype=int)
        return bool(0 < x.sum() < len(x))

Bases: CombinatorialProblem

Minimum-PMU full observability (zero-injection buses not modeled).

Variables: one placement bit per bus plus ceil(log2(deg_i + 1)) slack bits per bus. The IEEE 14-bus system needs 14 + 26 = 40 binaries — small for an annealer, instructive for counting why qubit budgets matter.

Source code in src/qugrid/problems/pmu.py

class PMUPlacement(CombinatorialProblem):
    """Minimum-PMU full observability (zero-injection buses not modeled).

    Variables: one placement bit per bus plus ``ceil(log2(deg_i + 1))`` slack
    bits per bus. The IEEE 14-bus system needs 14 + 26 = 40 binaries — small
    for an annealer, instructive for counting why qubit budgets matter.
    """

    def __init__(self, net: Network, penalty: float | None = None):
        self.net = net
        a = net.adjacency()
        self.neighbors = [np.flatnonzero(a[i]).tolist() for i in range(net.n_bus)]
        self.degrees = np.array([len(nb) for nb in self.neighbors])
        self.penalty = penalty if penalty is not None else float(2 * net.n_bus)

    @cached_property
    def qubo(self) -> QUBO:
        n = self.net.n_bus
        bld = QUBOBuilder()
        x = [bld.var(f"pmu@bus{int(self.net.bus[i, 0])}") for i in range(n)]
        for i in range(n):
            bld.add_linear(x[i], 1.0)  # objective: number of PMUs
        for i in range(n):
            # x_i + sum_{j in N(i)} x_j - 1 - slack_i = 0, slack_i in [0, deg_i]
            n_slack_bits = max(1, int(np.ceil(np.log2(self.degrees[i] + 1))))
            terms = [(x[i], 1.0)] + [(x[j], 1.0) for j in self.neighbors[i]]
            terms += [
                (bld.var(f"s[{i},{k}]"), -float(2**k)) for k in range(n_slack_bits)
            ]
            bld.add_squared_penalty(terms, -1.0, weight=self.penalty)
        return bld.build()

    def coverage(self, x: np.ndarray) -> np.ndarray:
        """Boolean per bus: observed by the placement bits in ``x``."""
        x = np.asarray(x, dtype=int)
        n = self.net.n_bus
        placed = x[:n].astype(bool)
        covered = placed.copy()
        for i in range(n):
            if placed[i]:
                covered[self.neighbors[i]] = True
        return covered

    def decode(self, x: np.ndarray) -> dict:
        n = self.net.n_bus
        placed = np.flatnonzero(np.asarray(x[:n], dtype=int))
        covered = self.coverage(x)
        return {
            "pmu_buses": [int(self.net.bus[i, 0]) for i in placed],
            "n_pmu": int(len(placed)),
            "n_uncovered": int((~covered).sum()),
            "covered": covered,
        }

    def is_feasible(self, x: np.ndarray) -> bool:
        return bool(self.coverage(x).all())

    def constraint_residual(self, x: np.ndarray) -> float:
        """Number of buses the placement bits in ``x`` leave unobserved."""
        return float((~self.coverage(x)).sum())

    def reference(self) -> dict:
        """Exact minimum dominating set by enumeration over placement bits.

        Slack bits do not need enumeration: for a feasible placement each
        coverage slack has one exact value (``coverage - 1``), so the full
        QUBO assignment — and with it a reference ``objective`` that
        :meth:`Result.gap` can use — is reconstructed in closed form.
        """
        n = self.net.n_bus
        if n > 20:
            raise ValueError("exact enumeration limited to 20 buses")
        best_x: np.ndarray | None = None
        best = n + 1
        for k in range(2**n):
            x = np.array([(k >> i) & 1 for i in range(n)], dtype=int)
            if x.sum() >= best:
                continue
            if self.coverage(x).all():
                best = int(x.sum())
                best_x = x
        assert best_x is not None
        bits = [best_x]
        for i in range(n):
            slack = int(best_x[i] + best_x[self.neighbors[i]].sum()) - 1
            n_slack_bits = max(1, int(np.ceil(np.log2(self.degrees[i] + 1))))
            bits.append(np.array([(slack >> k) & 1 for k in range(n_slack_bits)]))
        x_full = np.concatenate(bits).astype(int)
        out = self.decode(x_full)
        out["x"] = x_full
        out["objective"] = float(self.qubo.energy(x_full))
        return out

coverage(x)

Boolean per bus: observed by the placement bits in x.

Source code in src/qugrid/problems/pmu.py

def coverage(self, x: np.ndarray) -> np.ndarray:
    """Boolean per bus: observed by the placement bits in ``x``."""
    x = np.asarray(x, dtype=int)
    n = self.net.n_bus
    placed = x[:n].astype(bool)
    covered = placed.copy()
    for i in range(n):
        if placed[i]:
            covered[self.neighbors[i]] = True
    return covered

constraint_residual(x)

Number of buses the placement bits in x leave unobserved.

Source code in src/qugrid/problems/pmu.py

def constraint_residual(self, x: np.ndarray) -> float:
    """Number of buses the placement bits in ``x`` leave unobserved."""
    return float((~self.coverage(x)).sum())

reference()

Exact minimum dominating set by enumeration over placement bits.

Slack bits do not need enumeration: for a feasible placement each coverage slack has one exact value (coverage - 1), so the full QUBO assignment — and with it a reference objective that :meth:Result.gap can use — is reconstructed in closed form.

Source code in src/qugrid/problems/pmu.py

def reference(self) -> dict:
    """Exact minimum dominating set by enumeration over placement bits.

    Slack bits do not need enumeration: for a feasible placement each
    coverage slack has one exact value (``coverage - 1``), so the full
    QUBO assignment — and with it a reference ``objective`` that
    :meth:`Result.gap` can use — is reconstructed in closed form.
    """
    n = self.net.n_bus
    if n > 20:
        raise ValueError("exact enumeration limited to 20 buses")
    best_x: np.ndarray | None = None
    best = n + 1
    for k in range(2**n):
        x = np.array([(k >> i) & 1 for i in range(n)], dtype=int)
        if x.sum() >= best:
            continue
        if self.coverage(x).all():
            best = int(x.sum())
            best_x = x
    assert best_x is not None
    bits = [best_x]
    for i in range(n):
        slack = int(best_x[i] + best_x[self.neighbors[i]].sum()) - 1
        n_slack_bits = max(1, int(np.ceil(np.log2(self.degrees[i] + 1))))
        bits.append(np.array([(slack >> k) & 1 for k in range(n_slack_bits)]))
    x_full = np.concatenate(bits).astype(int)
    out = self.decode(x_full)
    out["x"] = x_full
    out["objective"] = float(self.qubo.energy(x_full))
    return out

The reduced DC power flow system B' theta = P as a solver-ready problem.

The slack bus row/column are removed (its angle is fixed), leaving a symmetric positive-definite system — exactly what HHL-style algorithms want. context carries everything needed to map x back to angles.

Source code in src/qugrid/problems/power_flow.py

def dc_power_flow(net: Network) -> LinearSystemProblem:
    """The reduced DC power flow system ``B' theta = P`` as a solver-ready problem.

    The slack bus row/column are removed (its angle is fixed), leaving a
    symmetric positive-definite system — exactly what HHL-style algorithms
    want. ``context`` carries everything needed to map ``x`` back to angles.
    """
    bbus, bf, pbusinj, pfinj = net.bdc()
    pbus = np.real(net.sbus()) - net.bus[:, idx.GS] / net.baseMVA
    ref = net.ref
    keep = np.setdiff1d(np.arange(net.n_bus), [ref])
    theta_ref = float(np.deg2rad(net.bus[ref, idx.VA]))
    a = bbus[np.ix_(keep, keep)]
    b = pbus[keep] - pbusinj[keep] - bbus[keep, ref] * theta_ref
    names = [f"theta_bus{int(net.bus[i, 0])}" for i in keep]
    return LinearSystemProblem(
        a=a,
        b=b,
        names=names,
        unit="rad",
        context={"net": net, "ref": ref, "keep": keep, "theta_ref": theta_ref, "bf": bf,
                 "pfinj": pfinj},
    )

Newton-Raphson AC power flow with a caller-supplied linear solver.

solve_fn(problem: LinearSystemProblem) -> np.ndarray receives each Newton step J dx = -F and returns dx. Pass an exact solver, HHL, or VQLS — this is the standard hybrid quantum-classical decomposition of AC power flow.

Source code in src/qugrid/problems/power_flow.py

def newton_with_linear_solver(
    net: Network,
    solve_fn,
    tol: float = 1e-6,
    max_iter: int = 20,
) -> HybridNewtonResult:
    """Newton-Raphson AC power flow with a caller-supplied linear solver.

    ``solve_fn(problem: LinearSystemProblem) -> np.ndarray`` receives each
    Newton step ``J dx = -F`` and returns ``dx``. Pass an exact solver, HHL,
    or VQLS — this is the standard hybrid quantum-classical decomposition of
    AC power flow.
    """
    ybus = net.ybus()
    sbus = net.sbus()
    ref, pv, pq = net.ref, net.pv, net.pq
    del ref
    pvpq = np.concatenate([pv, pq]).astype(int)

    vm = net.bus[:, idx.VM].copy()
    va = np.deg2rad(net.bus[:, idx.VA])
    on = net.gen_on
    vm[net.gen_bus[on]] = net.gen[on, idx.VG]
    v = vm * np.exp(1j * va)

    history: list[float] = []
    inner_residuals: list[float] = []
    converged = False
    it = 0
    for it in range(1, max_iter + 1):
        ibus = ybus @ v
        mis = v * np.conj(ibus) - sbus
        f_vec = np.concatenate([np.real(mis[pvpq]), np.imag(mis[pq])])
        norm = float(np.linalg.norm(f_vec, np.inf)) if f_vec.size else 0.0
        history.append(norm)
        if norm < tol:
            converged = True
            break

        diag_v = np.diag(v)
        diag_i = np.diag(ibus)
        diag_vn = np.diag(v / np.abs(v))
        ds_dva = 1j * diag_v @ np.conj(diag_i - ybus @ diag_v)
        ds_dvm = diag_v @ np.conj(ybus @ diag_vn) + np.conj(diag_i) @ diag_vn
        j11 = np.real(ds_dva[np.ix_(pvpq, pvpq)])
        j12 = np.real(ds_dvm[np.ix_(pvpq, pq)])
        j21 = np.imag(ds_dva[np.ix_(pq, pvpq)])
        j22 = np.imag(ds_dvm[np.ix_(pq, pq)])
        jac = np.block([[j11, j12], [j21, j22]])

        step = LinearSystemProblem(
            a=jac, b=-f_vec, unit="mixed rad/pu", context={"iteration": it}
        )
        dx = np.asarray(solve_fn(step), dtype=float)
        inner_residuals.append(float(np.linalg.norm(jac @ dx + f_vec)))

        n1 = len(pvpq)
        va[pvpq] += dx[:n1]
        vm[pq] += dx[n1:]
        v = vm * np.exp(1j * va)

    s_from, s_to = _branch_flows(net, v)
    ac = ACResult(
        vm=np.abs(v),
        va=np.angle(v),
        converged=converged,
        iterations=it,
        mismatch_history=history,
        s_from_mva=s_from,
        s_to_mva=s_to,
    )
    return HybridNewtonResult(ac=ac, linear_solves=len(inner_residuals),
                              inner_residuals=inner_residuals)

Sample load patterns and label them secure/insecure by DC power flow.

Each load bus gets an independent multiplier drawn uniformly from load_range; generation is scaled proportionally to keep balance. The label is 1 when the maximum branch loading (|flow| / RATE_A) exceeds threshold in the base case or in any screened contingency.

contingencies: "n-1" (default) screens every single-branch outage that keeps the grid connected; None labels on the base case only; a list of branch indices screens exactly those outages.

Features are the multipliers — low-dimensional (one per load bus) and physically meaningful, sized for today's quantum feature maps.

Source code in src/qugrid/problems/screening.py

def screening_dataset(
    net: Network,
    n_samples: int = 200,
    load_range: tuple[float, float] = (0.6, 1.4),
    threshold: float = 1.0,
    contingencies: list[int] | str | None = "n-1",
    seed: int = 0,
) -> ScreeningDataset:
    """Sample load patterns and label them secure/insecure by DC power flow.

    Each load bus gets an independent multiplier drawn uniformly from
    ``load_range``; generation is scaled proportionally to keep balance.
    The label is 1 when the maximum branch loading (|flow| / RATE_A) exceeds
    ``threshold`` in the base case **or in any screened contingency**.

    ``contingencies``: ``"n-1"`` (default) screens every single-branch outage
    that keeps the grid connected; ``None`` labels on the base case only; a
    list of branch indices screens exactly those outages.

    Features are the multipliers — low-dimensional (one per load bus) and
    physically meaningful, sized for today's quantum feature maps.
    """
    if not np.any(net.branch[:, idx.RATE_A] > 0):
        raise ValueError("network has no branch ratings (RATE_A); cannot label security")
    if contingencies == "n-1":
        cont = screenable_contingencies(net)
    elif contingencies is None:
        cont = []
    else:
        cont = list(contingencies)
        for k in cont:
            if not _connected_without(net, k):
                raise ValueError(f"outage of branch {k} disconnects the network")
    topologies = [net.copy()] + [net.drop_branch(k) for k in cont]

    load_buses = np.flatnonzero(net.bus[:, idx.PD] > 0)
    rng = np.random.default_rng(seed)
    x = rng.uniform(load_range[0], load_range[1], size=(n_samples, len(load_buses)))
    y = np.zeros(n_samples, dtype=int)
    base_load = net.bus[:, idx.PD].sum()
    for s in range(n_samples):
        worst = 0.0
        for topo in topologies:
            sample = topo.copy()
            factors = np.ones(sample.n_bus)
            factors[load_buses] = x[s]
            sample.bus[:, idx.PD] *= factors
            sample.bus[:, idx.QD] *= factors
            # keep generation-load balance by scaling PG proportionally
            scale = sample.bus[:, idx.PD].sum() / max(base_load, 1e-9)
            sample.gen[:, idx.PG] *= scale
            worst = max(worst, solve_dc(sample).max_loading())
            if worst > threshold:
                break
        y[s] = int(worst > threshold)
    names = [f"load_bus{int(net.bus[i, idx.BUS_I])}" for i in load_buses]
    return ScreeningDataset(
        x=x,
        y=y,
        feature_names=names,
        threshold=threshold,
        contingencies=cont,
        net_name=net.name,
    )

Encodings

Quadratic unconstrained binary optimization: minimize x^T Q x + offset.

Source code in src/qugrid/problems/base.py

@dataclass
class QUBO:
    """Quadratic unconstrained binary optimization: minimize ``x^T Q x + offset``."""

    q: np.ndarray  # (n, n); symmetrized on construction
    offset: float = 0.0
    names: list[str] = field(default_factory=list)

    def __post_init__(self) -> None:
        q = np.asarray(self.q, dtype=float)
        self.q = 0.5 * (q + q.T)
        if not self.names:
            self.names = [f"x{i}" for i in range(self.n)]

    @property
    def n(self) -> int:
        return self.q.shape[0]

    def energy(self, x: np.ndarray) -> float | np.ndarray:
        """Objective value of bitstring(s); ``x`` is (n,) or (batch, n) of 0/1."""
        x = np.asarray(x, dtype=float)
        if x.ndim == 1:
            return float(x @ self.q @ x + self.offset)
        return np.einsum("bi,ij,bj->b", x, self.q, x) + self.offset

    def to_ising(self) -> Ising:
        """Exact change of variables ``s = 1 - 2x``. Energies are preserved."""
        d = np.diag(self.q).copy()
        off = self.q - np.diag(d)  # symmetric, zero diagonal
        row = off.sum(axis=1)
        j = np.triu(off / 2.0, k=1)  # pair (i<j) contributes 2*Q_ij*x_i*x_j -> J_ij = Q_ij/2
        h = -row / 2.0 - d / 2.0
        offset = self.offset + off.sum() / 4.0 + d.sum() / 2.0
        return Ising(h=h, j=j, offset=offset)

    def bits_from_index(self, k: int) -> np.ndarray:
        return np.array([(k >> i) & 1 for i in range(self.n)], dtype=int)

    def dynamic_range(self, db: bool = False) -> float:
        """``max|coef| / min nonzero |coef|`` over the energy coefficients.

        Coefficients are the linear terms ``Q_ii`` and the pair couplings
        ``2 Q_ij`` (i < j); the offset shifts every energy equally and is
        excluded. Penalty-folded formulations stretch this ratio, and
        sampling solvers stop resolving the costs underneath: the 2-unit
        unit commitment study's penalized QUBO measures 3.9e3 here, and
        example 07 documents the QAOA success-probability collapse on it.
        ``db=True`` returns ``20 log10`` of the ratio. A QUBO with no
        nonzero coefficient reports 1 (0 dB).
        """
        lin = np.abs(np.diag(self.q))
        pairs = 2.0 * np.abs(self.q[np.triu_indices(self.n, k=1)])
        coefs = np.concatenate([lin, pairs])
        nonzero = coefs[coefs > 0.0]
        if nonzero.size == 0:
            return 0.0 if db else 1.0
        ratio = float(nonzero.max() / nonzero.min())
        return float(20.0 * np.log10(ratio)) if db else ratio

energy(x)

Objective value of bitstring(s); x is (n,) or (batch, n) of 0/1.

Source code in src/qugrid/problems/base.py

def energy(self, x: np.ndarray) -> float | np.ndarray:
    """Objective value of bitstring(s); ``x`` is (n,) or (batch, n) of 0/1."""
    x = np.asarray(x, dtype=float)
    if x.ndim == 1:
        return float(x @ self.q @ x + self.offset)
    return np.einsum("bi,ij,bj->b", x, self.q, x) + self.offset

to_ising()

Exact change of variables s = 1 - 2x. Energies are preserved.

Source code in src/qugrid/problems/base.py

def to_ising(self) -> Ising:
    """Exact change of variables ``s = 1 - 2x``. Energies are preserved."""
    d = np.diag(self.q).copy()
    off = self.q - np.diag(d)  # symmetric, zero diagonal
    row = off.sum(axis=1)
    j = np.triu(off / 2.0, k=1)  # pair (i<j) contributes 2*Q_ij*x_i*x_j -> J_ij = Q_ij/2
    h = -row / 2.0 - d / 2.0
    offset = self.offset + off.sum() / 4.0 + d.sum() / 2.0
    return Ising(h=h, j=j, offset=offset)

dynamic_range(db=False)

max|coef| / min nonzero |coef| over the energy coefficients.

Coefficients are the linear terms Q_ii and the pair couplings 2 Q_ij (i < j); the offset shifts every energy equally and is excluded. Penalty-folded formulations stretch this ratio, and sampling solvers stop resolving the costs underneath: the 2-unit unit commitment study's penalized QUBO measures 3.9e3 here, and example 07 documents the QAOA success-probability collapse on it. db=True returns 20 log10 of the ratio. A QUBO with no nonzero coefficient reports 1 (0 dB).

Source code in src/qugrid/problems/base.py

def dynamic_range(self, db: bool = False) -> float:
    """``max|coef| / min nonzero |coef|`` over the energy coefficients.

    Coefficients are the linear terms ``Q_ii`` and the pair couplings
    ``2 Q_ij`` (i < j); the offset shifts every energy equally and is
    excluded. Penalty-folded formulations stretch this ratio, and
    sampling solvers stop resolving the costs underneath: the 2-unit
    unit commitment study's penalized QUBO measures 3.9e3 here, and
    example 07 documents the QAOA success-probability collapse on it.
    ``db=True`` returns ``20 log10`` of the ratio. A QUBO with no
    nonzero coefficient reports 1 (0 dB).
    """
    lin = np.abs(np.diag(self.q))
    pairs = 2.0 * np.abs(self.q[np.triu_indices(self.n, k=1)])
    coefs = np.concatenate([lin, pairs])
    nonzero = coefs[coefs > 0.0]
    if nonzero.size == 0:
        return 0.0 if db else 1.0
    ratio = float(nonzero.max() / nonzero.min())
    return float(20.0 * np.log10(ratio)) if db else ratio

Ising Hamiltonian H(s) = sum_{i<j} J_ij s_i s_j + sum_i h_i s_i + offset.

Source code in src/qugrid/problems/base.py

@dataclass
class Ising:
    """Ising Hamiltonian ``H(s) = sum_{i<j} J_ij s_i s_j + sum_i h_i s_i + offset``."""

    h: np.ndarray  # (n,)
    j: np.ndarray  # (n, n), strictly upper triangular (J_ij defined for i < j)
    offset: float = 0.0

    def __post_init__(self) -> None:
        # Normalize once so every consumer (SA local fields, adapters, QAOA)
        # can rely on the strictly-upper-triangular invariant.
        self.j = np.triu(np.asarray(self.j, dtype=float), k=1)
        self.h = np.asarray(self.h, dtype=float)

    @property
    def n(self) -> int:
        return len(self.h)

    def energy(self, s: np.ndarray) -> float | np.ndarray:
        """Energy of spin configuration(s); ``s`` is (n,) or (batch, n) of +-1."""
        s = np.asarray(s, dtype=float)
        if s.ndim == 1:
            return float(s @ self.j @ s + self.h @ s + self.offset)
        return np.einsum("bi,ij,bj->b", s, self.j, s) + s @ self.h + self.offset

    def all_energies(self) -> np.ndarray:
        """Energies of all ``2^n`` configurations (index bit ``i`` = variable ``i``).

        Memory stays at O(2^n) — spins are generated per variable, not cached
        per variable, so n = 24 costs ~400 MB peak instead of ~3 GB.
        """
        if self.n > 24:
            raise ValueError(f"refusing to enumerate 2^{self.n} states")
        states = np.arange(2**self.n, dtype=np.int64)

        def spin(i: int) -> np.ndarray:
            return 1.0 - 2.0 * ((states >> i) & 1)  # bit 0 -> spin +1

        energies = np.full(2**self.n, self.offset, dtype=float)
        for i in range(self.n):
            if self.h[i] != 0.0:
                energies += self.h[i] * spin(i)
        ii, jj = np.nonzero(self.j)
        last_a, s_a = -1, None
        for a, b in zip(ii, jj):
            if a != last_a:  # ii is row-sorted; cache the left spin per row
                s_a, last_a = spin(a), a
            energies += (self.j[a, b] * s_a) * spin(b)
        return energies

energy(s)

Energy of spin configuration(s); s is (n,) or (batch, n) of +-1.

Source code in src/qugrid/problems/base.py

def energy(self, s: np.ndarray) -> float | np.ndarray:
    """Energy of spin configuration(s); ``s`` is (n,) or (batch, n) of +-1."""
    s = np.asarray(s, dtype=float)
    if s.ndim == 1:
        return float(s @ self.j @ s + self.h @ s + self.offset)
    return np.einsum("bi,ij,bj->b", s, self.j, s) + s @ self.h + self.offset

all_energies()

Energies of all 2^n configurations (index bit i = variable i).

Memory stays at O(2^n) — spins are generated per variable, not cached per variable, so n = 24 costs ~400 MB peak instead of ~3 GB.

Source code in src/qugrid/problems/base.py

def all_energies(self) -> np.ndarray:
    """Energies of all ``2^n`` configurations (index bit ``i`` = variable ``i``).

    Memory stays at O(2^n) — spins are generated per variable, not cached
    per variable, so n = 24 costs ~400 MB peak instead of ~3 GB.
    """
    if self.n > 24:
        raise ValueError(f"refusing to enumerate 2^{self.n} states")
    states = np.arange(2**self.n, dtype=np.int64)

    def spin(i: int) -> np.ndarray:
        return 1.0 - 2.0 * ((states >> i) & 1)  # bit 0 -> spin +1

    energies = np.full(2**self.n, self.offset, dtype=float)
    for i in range(self.n):
        if self.h[i] != 0.0:
            energies += self.h[i] * spin(i)
    ii, jj = np.nonzero(self.j)
    last_a, s_a = -1, None
    for a, b in zip(ii, jj):
        if a != last_a:  # ii is row-sorted; cache the left spin per row
            s_a, last_a = spin(a), a
        energies += (self.j[a, b] * s_a) * spin(b)
    return energies

A linear system A x = b destined for a (quantum) linear solver.

scale and names carry the engineering meaning of x so that solvers can hand back answers in physical units.

Source code in src/qugrid/problems/base.py

@dataclass
class LinearSystemProblem:
    """A linear system ``A x = b`` destined for a (quantum) linear solver.

    ``scale`` and ``names`` carry the engineering meaning of ``x`` so that
    solvers can hand back answers in physical units.
    """

    a: np.ndarray
    b: np.ndarray
    names: list[str] = field(default_factory=list)
    unit: str = ""
    context: dict = field(default_factory=dict)

    def __post_init__(self) -> None:
        self.a = np.asarray(self.a, dtype=float)
        self.b = np.asarray(self.b, dtype=float)
        if not self.names:
            self.names = [f"x{i}" for i in range(len(self.b))]

    @property
    def n(self) -> int:
        return len(self.b)

    @property
    def is_hermitian(self) -> bool:
        return bool(np.allclose(self.a, self.a.T))

    def condition_number(self) -> float:
        return float(np.linalg.cond(self.a))

    def solve_exact(self) -> np.ndarray:
        return np.linalg.solve(self.a, self.b)

    def padded(self) -> tuple[np.ndarray, np.ndarray, int]:
        """Embed into the next power-of-two dimension (identity padding).

        Quantum registers hold ``2^m`` amplitudes; the padding is inert
        (``A'`` block diagonal with the identity, ``b'`` zero on the pad).
        Returns ``(A', b', original_dimension)``.
        """
        n = self.n
        m = 1 if n <= 1 else int(np.ceil(np.log2(n)))
        full = 2**m
        a = np.eye(full)
        a[:n, :n] = self.a
        b = np.zeros(full)
        b[:n] = self.b
        return a, b, n

    def __repr__(self) -> str:
        return (
            f"LinearSystemProblem(n={self.n}, unit={self.unit!r}, "
            f"cond={self.condition_number():.3g})"
        )

padded()

Embed into the next power-of-two dimension (identity padding).

Quantum registers hold 2^m amplitudes; the padding is inert (A' block diagonal with the identity, b' zero on the pad). Returns (A', b', original_dimension).

Source code in src/qugrid/problems/base.py

def padded(self) -> tuple[np.ndarray, np.ndarray, int]:
    """Embed into the next power-of-two dimension (identity padding).

    Quantum registers hold ``2^m`` amplitudes; the padding is inert
    (``A'`` block diagonal with the identity, ``b'`` zero on the pad).
    Returns ``(A', b', original_dimension)``.
    """
    n = self.n
    m = 1 if n <= 1 else int(np.ceil(np.log2(n)))
    full = 2**m
    a = np.eye(full)
    a[:n, :n] = self.a
    b = np.zeros(full)
    b[:n] = self.b
    return a, b, n

Source code in src/qugrid/problems/builder.py

class QUBOBuilder:
    def __init__(self) -> None:
        self._names: list[str] = []
        self._index: dict[str, int] = {}
        self._lin: dict[int, float] = {}
        self._quad: dict[tuple[int, int], float] = {}
        self._const: float = 0.0
        self._n_ineq: int = 0

    def var(self, name: str) -> int:
        """Register a binary variable and return its index (idempotent)."""
        if name not in self._index:
            self._index[name] = len(self._names)
            self._names.append(name)
        return self._index[name]

    @property
    def n(self) -> int:
        return len(self._names)

    def add_constant(self, c: float) -> None:
        self._const += float(c)

    def add_linear(self, i: int, c: float) -> None:
        self._lin[i] = self._lin.get(i, 0.0) + float(c)

    def add_quadratic(self, i: int, j: int, c: float) -> None:
        if i == j:
            # x_i^2 == x_i for binaries
            self.add_linear(i, c)
            return
        key = (min(i, j), max(i, j))
        self._quad[key] = self._quad.get(key, 0.0) + float(c)

    def add_squared_penalty(
        self, terms: list[tuple[int, float]], constant: float, weight: float
    ) -> None:
        """Add ``weight * (sum coef_i x_i + constant)^2``, expanded exactly."""
        w = float(weight)
        c = float(constant)
        self.add_constant(w * c * c)
        for i, a in terms:
            self.add_linear(i, w * (a * a + 2.0 * a * c))
        for k in range(len(terms)):
            i, a = terms[k]
            for m in range(k + 1, len(terms)):
                j, b = terms[m]
                self.add_quadratic(i, j, 2.0 * w * a * b)

    def add_inequality(
        self,
        terms: list[tuple[int, float]],
        constant: float,
        weight: float = 1.0,
        method: str = "slack",
        lam: tuple[float, float] | None = None,
        name: str | None = None,
    ) -> None:
        """Penalize the constraint ``sum_i coef_i * x_i + constant >= 0``.

        ``method="slack"`` (default) is the textbook route (Lucas,
        arXiv:1302.5843): register binary slack bits covering ``[0, U]``
        with ``U`` the maximum of the left side over all bitstrings, and add
        ``weight * (sum_i coef_i x_i + constant - S)^2``. Exact whenever the
        left side is integer-valued; costs ``ceil(log2(U + 1))`` extra bits
        per constraint.

        ``method="unbalanced"`` is unbalanced penalization
        (Montanez-Barrera et al., Quantum Sci. Technol. 9, 025022, 2024;
        arXiv:2211.13914): add ``-l1*g + l2*g^2`` for
        ``g = sum_i coef_i x_i + constant``, with ``lam=(l1, l2)``. Zero
        slack bits, so the qubit count stops growing with the number of
        inequalities — but the encoding is a heuristic: the unconstrained
        minimum coincides with the constrained optimum only for suitable
        ``(l1, l2)``, which need tuning per problem family. ``weight`` is
        ignored on this path.
        """
        if name is None:
            name = f"ineq{self._n_ineq}"
        self._n_ineq += 1
        if method == "slack":
            g_max = float(constant) + sum(max(a, 0.0) for _, a in terms)
            u = int(round(g_max))
            if abs(g_max - u) > 1e-9:
                raise ValueError("slack encoding needs an integer-valued left side")
            if u < 0:
                raise ValueError("constraint can never hold: left side is at most negative")
            n_bits = max(1, int(np.ceil(np.log2(u + 1)))) if u > 0 else 0
            slack = [(self.var(f"{name}_s{k}"), -float(2**k)) for k in range(n_bits)]
            self.add_squared_penalty(list(terms) + slack, constant, weight)
        elif method == "unbalanced":
            l1, l2 = lam if lam is not None else (1.0, 1.0)
            self.add_constant(-l1 * float(constant))
            for i, a in terms:
                self.add_linear(i, -l1 * a)
            self.add_squared_penalty(list(terms), constant, weight=l2)
        else:
            raise ValueError(f"unknown method {method!r}; use 'slack' or 'unbalanced'")

    def build(self) -> QUBO:
        q = np.zeros((self.n, self.n))
        for i, c in self._lin.items():
            q[i, i] += c
        for (i, j), c in self._quad.items():
            q[i, j] += c / 2.0
            q[j, i] += c / 2.0
        return QUBO(q=q, offset=self._const, names=list(self._names))

var(name)

Register a binary variable and return its index (idempotent).

Source code in src/qugrid/problems/builder.py

def var(self, name: str) -> int:
    """Register a binary variable and return its index (idempotent)."""
    if name not in self._index:
        self._index[name] = len(self._names)
        self._names.append(name)
    return self._index[name]

add_squared_penalty(terms, constant, weight)

Add weight * (sum coef_i x_i + constant)^2, expanded exactly.

Source code in src/qugrid/problems/builder.py

def add_squared_penalty(
    self, terms: list[tuple[int, float]], constant: float, weight: float
) -> None:
    """Add ``weight * (sum coef_i x_i + constant)^2``, expanded exactly."""
    w = float(weight)
    c = float(constant)
    self.add_constant(w * c * c)
    for i, a in terms:
        self.add_linear(i, w * (a * a + 2.0 * a * c))
    for k in range(len(terms)):
        i, a = terms[k]
        for m in range(k + 1, len(terms)):
            j, b = terms[m]
            self.add_quadratic(i, j, 2.0 * w * a * b)

add_inequality(terms, constant, weight=1.0, method='slack', lam=None, name=None)

Penalize the constraint sum_i coef_i * x_i + constant >= 0.

method="slack" (default) is the textbook route (Lucas, arXiv:1302.5843): register binary slack bits covering [0, U] with U the maximum of the left side over all bitstrings, and add weight * (sum_i coef_i x_i + constant - S)^2. Exact whenever the left side is integer-valued; costs ceil(log2(U + 1)) extra bits per constraint.

method="unbalanced" is unbalanced penalization (Montanez-Barrera et al., Quantum Sci. Technol. 9, 025022, 2024; arXiv:2211.13914): add -l1*g + l2*g^2 for g = sum_i coef_i x_i + constant, with lam=(l1, l2). Zero slack bits, so the qubit count stops growing with the number of inequalities — but the encoding is a heuristic: the unconstrained minimum coincides with the constrained optimum only for suitable (l1, l2), which need tuning per problem family. weight is ignored on this path.

Source code in src/qugrid/problems/builder.py

def add_inequality(
    self,
    terms: list[tuple[int, float]],
    constant: float,
    weight: float = 1.0,
    method: str = "slack",
    lam: tuple[float, float] | None = None,
    name: str | None = None,
) -> None:
    """Penalize the constraint ``sum_i coef_i * x_i + constant >= 0``.

    ``method="slack"`` (default) is the textbook route (Lucas,
    arXiv:1302.5843): register binary slack bits covering ``[0, U]``
    with ``U`` the maximum of the left side over all bitstrings, and add
    ``weight * (sum_i coef_i x_i + constant - S)^2``. Exact whenever the
    left side is integer-valued; costs ``ceil(log2(U + 1))`` extra bits
    per constraint.

    ``method="unbalanced"`` is unbalanced penalization
    (Montanez-Barrera et al., Quantum Sci. Technol. 9, 025022, 2024;
    arXiv:2211.13914): add ``-l1*g + l2*g^2`` for
    ``g = sum_i coef_i x_i + constant``, with ``lam=(l1, l2)``. Zero
    slack bits, so the qubit count stops growing with the number of
    inequalities — but the encoding is a heuristic: the unconstrained
    minimum coincides with the constrained optimum only for suitable
    ``(l1, l2)``, which need tuning per problem family. ``weight`` is
    ignored on this path.
    """
    if name is None:
        name = f"ineq{self._n_ineq}"
    self._n_ineq += 1
    if method == "slack":
        g_max = float(constant) + sum(max(a, 0.0) for _, a in terms)
        u = int(round(g_max))
        if abs(g_max - u) > 1e-9:
            raise ValueError("slack encoding needs an integer-valued left side")
        if u < 0:
            raise ValueError("constraint can never hold: left side is at most negative")
        n_bits = max(1, int(np.ceil(np.log2(u + 1)))) if u > 0 else 0
        slack = [(self.var(f"{name}_s{k}"), -float(2**k)) for k in range(n_bits)]
        self.add_squared_penalty(list(terms) + slack, constant, weight)
    elif method == "unbalanced":
        l1, l2 = lam if lam is not None else (1.0, 1.0)
        self.add_constant(-l1 * float(constant))
        for i, a in terms:
            self.add_linear(i, -l1 * a)
        self.add_squared_penalty(list(terms), constant, weight=l2)
    else:
        raise ValueError(f"unknown method {method!r}; use 'slack' or 'unbalanced'")

Results and benchmarking

Outcome of solving a QuGrid problem.

Source code in src/qugrid/solvers/base.py

@dataclass
class Result:
    """Outcome of solving a QuGrid problem."""

    solver: str
    problem: Any
    x: np.ndarray | None = None
    objective: float | None = None
    decoded: dict = field(default_factory=dict)
    feasible: bool | None = None
    history: list[float] = field(default_factory=list)
    top_states: list[tuple[str, float, float]] = field(default_factory=list)
    # (bitstring, probability, energy) for combinatorial quantum solvers
    resources: dict = field(default_factory=dict)
    reference: dict | None = None
    extras: dict = field(default_factory=dict)

    # ------------------------------------------------------------ derived
    def gap(self) -> float | None:
        """Relative optimality gap versus the reference objective (if known)."""
        if self.reference is None or self.objective is None:
            return None
        ref = self.reference.get("objective")
        if ref is None:
            return None
        denom = max(abs(ref), 1e-12)
        return float((self.objective - ref) / denom)

    def success_probability(self) -> float | None:
        """Probability mass on the reference optimum (quantum solvers only)."""
        if self.reference is None or not self.top_states:
            return None
        ref = self.reference.get("objective")
        if ref is None:
            return None
        tol = 1e-9 * max(1.0, abs(ref))
        return float(
            sum(p for _, p, e in self.top_states if abs(e - ref) <= tol)
        )

    # ------------------------------------------------------------- display
    def summary(self) -> str:
        lines = [f"QuGrid result | solver={self.solver} | {self._problem_name()}"]
        if self.objective is not None:
            lines.append(f"  objective          {self.objective:,.6g}")
        if self.feasible is not None:
            lines.append(f"  feasible           {'yes' if self.feasible else 'NO'}")
        g = self.gap()
        if g is not None:
            if abs(g) < 1e-9:  # display float noise as the exact match it is
                g = 0.0
            lines.append(f"  gap vs reference   {100 * g:.3g}%")
        sp = self.success_probability()
        if sp is not None:
            lines.append(f"  P(optimum)         {sp:.3f}")
        pr = self.extras.get("p_optimum_repaired")
        if pr is not None:
            lines.append(f"  P(opt | repaired)  {pr:.3f}")
        for key, val in self.decoded.items():
            if isinstance(val, (int, float, np.floating)):
                lines.append(f"  {key:<18} {val:,.6g}")
            elif isinstance(val, tuple) and all(isinstance(v, (int, float)) for v in val):
                shown = tuple(v if isinstance(v, (bool, int)) else round(float(v), 4) for v in val)
                lines.append(f"  {key:<18} {shown}")
        if self.resources:
            res = ", ".join(f"{k}={v}" for k, v in self.resources.items())
            lines.append(f"  resources          {res}")
        return "\n".join(lines)

    def _problem_name(self) -> str:
        p = self.problem
        name = type(p).__name__
        n = getattr(p, "n", None)
        return f"{name}(n={n})" if n is not None else name

    def plot_convergence(self, ax=None):
        """Optimizer trajectory (thin wrapper over :func:`qugrid.viz.plot_convergence`)."""
        from qugrid.viz import plot_convergence

        return plot_convergence(self, ax=ax)

    def __repr__(self) -> str:
        obj = f"{self.objective:.6g}" if self.objective is not None else "n/a"
        return f"Result(solver={self.solver!r}, objective={obj}, feasible={self.feasible})"

gap()

Relative optimality gap versus the reference objective (if known).

Source code in src/qugrid/solvers/base.py

def gap(self) -> float | None:
    """Relative optimality gap versus the reference objective (if known)."""
    if self.reference is None or self.objective is None:
        return None
    ref = self.reference.get("objective")
    if ref is None:
        return None
    denom = max(abs(ref), 1e-12)
    return float((self.objective - ref) / denom)

success_probability()

Probability mass on the reference optimum (quantum solvers only).

Source code in src/qugrid/solvers/base.py

def success_probability(self) -> float | None:
    """Probability mass on the reference optimum (quantum solvers only)."""
    if self.reference is None or not self.top_states:
        return None
    ref = self.reference.get("objective")
    if ref is None:
        return None
    tol = 1e-9 * max(1.0, abs(ref))
    return float(
        sum(p for _, p, e in self.top_states if abs(e - ref) <= tol)
    )

plot_convergence(ax=None)

Optimizer trajectory (thin wrapper over :func:qugrid.viz.plot_convergence).

Source code in src/qugrid/solvers/base.py

def plot_convergence(self, ax=None):
    """Optimizer trajectory (thin wrapper over :func:`qugrid.viz.plot_convergence`)."""
    from qugrid.viz import plot_convergence

    return plot_convergence(self, ax=ax)

Paper-ready benchmarking: sweep, aggregate, export.

The loop every quantum-for-power paper runs — problems x solvers x seeds — in one call, with the outputs reviewers ask for: a tidy DataFrame, a LaTeX table, and a JSON config that pins every seed and parameter.

import qugrid as qg from qugrid import bench problems = {"case9": qg.problems.Islanding(qg.cases.case9())} df = bench.sweep(problems, solvers=["exact", "sa", "qaoa"], seeds=range(3)) print(bench.summarize(df)) bench.save_run(df, "runs/islanding") # results.csv + summary.tex + config.json

sweep(problems, solvers, seeds=range(3), verbose=True, **common_kwargs)

Run every solver on every problem with every seed.

solvers is a list of registry names, or a dict mapping name to extra keyword arguments (for example {"qaoa": {"p": 3}}). Deterministic solvers still run once per seed; their spread doubles as a sanity check.

Source code in src/qugrid/bench.py

def sweep(
    problems: dict[str, object],
    solvers: list[str] | dict[str, dict],
    seeds=range(3),
    verbose: bool = True,
    **common_kwargs,
) -> pd.DataFrame:
    """Run every solver on every problem with every seed.

    ``solvers`` is a list of registry names, or a dict mapping name to extra
    keyword arguments (for example ``{"qaoa": {"p": 3}}``). Deterministic
    solvers still run once per seed; their spread doubles as a sanity check.
    """
    solver_items = (
        [(name, {}) for name in solvers]
        if isinstance(solvers, list)
        else list(solvers.items())
    )
    rows = []
    for label, problem in problems.items():
        for solver_name, extra in solver_items:
            for seed in seeds:
                kwargs = {**common_kwargs, **extra}
                res = solve(problem, solver=solver_name, seed=int(seed), **kwargs)
                gap = res.gap()
                row = {
                    "problem": label,
                    "solver": solver_name,
                    "seed": int(seed),
                    "objective": res.objective,
                    "gap": gap,
                    "feasible": res.feasible,
                    "success_probability": res.success_probability(),
                    "wall_time_s": res.resources.get("wall_time_s"),
                    "n_qubits": res.resources.get("n_qubits"),
                    "evaluations": res.resources.get("evaluations"),
                }
                rows.append(row)
                if verbose:
                    gap_s = f"{100 * gap:.2f}%" if gap is not None else "n/a"
                    print(
                        f"[bench] {label:<12} {solver_name:<7} seed={seed} "
                        f"objective={res.objective:,.4g} gap={gap_s}"
                    )
    return pd.DataFrame(rows)

summarize(df)

Mean and spread per (problem, solver), the table most papers print.

Source code in src/qugrid/bench.py

def summarize(df: pd.DataFrame) -> pd.DataFrame:
    """Mean and spread per (problem, solver), the table most papers print."""
    agg = df.groupby(["problem", "solver"]).agg(
        objective_mean=("objective", "mean"),
        objective_std=("objective", "std"),
        gap_mean=("gap", "mean"),
        feasible_rate=("feasible", "mean"),
        p_success_mean=("success_probability", "mean"),
        time_mean_s=("wall_time_s", "mean"),
        runs=("seed", "count"),
    )
    return agg.reset_index()

to_latex(df, caption='QuGrid benchmark', label='tab:qugrid')

Booktabs LaTeX for a summarized DataFrame (drop it straight into a paper).

Written by hand rather than through DataFrame.to_latex so the core install needs no jinja2 (pandas >= 2.0 routes to_latex through its Styler, which imports jinja2 at call time).

Source code in src/qugrid/bench.py

def to_latex(df: pd.DataFrame, caption: str = "QuGrid benchmark", label: str = "tab:qugrid"):
    """Booktabs LaTeX for a summarized DataFrame (drop it straight into a paper).

    Written by hand rather than through ``DataFrame.to_latex`` so the core
    install needs no jinja2 (pandas >= 2.0 routes ``to_latex`` through its
    Styler, which imports jinja2 at call time).
    """
    summary = df if "objective_mean" in df.columns else summarize(df)
    cols = list(summary.columns)
    align = "".join("l" if summary[c].dtype == object else "r" for c in cols)
    lines = [
        "\\begin{table}[t]",
        "\\centering",
        f"\\caption{{{caption}}}",
        f"\\label{{{label}}}",
        f"\\begin{{tabular}}{{{align}}}",
        "\\toprule",
        " & ".join(_latex_cell(c) for c in cols) + " \\\\",
        "\\midrule",
    ]
    for _, row in summary.iterrows():
        lines.append(" & ".join(_latex_cell(row[c]) for c in cols) + " \\\\")
    lines += ["\\bottomrule", "\\end{tabular}", "\\end{table}", ""]
    return "\n".join(lines)

save_run(df, out_dir, config=None, caption='QuGrid benchmark')

Write results.csv, summary.csv, summary.tex, config.json.

The config records package version, platform, and timestamp so the run can be cited and reproduced.

Source code in src/qugrid/bench.py

def save_run(
    df: pd.DataFrame,
    out_dir: str | Path,
    config: dict | None = None,
    caption: str = "QuGrid benchmark",
) -> Path:
    """Write ``results.csv``, ``summary.csv``, ``summary.tex``, ``config.json``.

    The config records package version, platform, and timestamp so the run
    can be cited and reproduced.
    """
    out = Path(out_dir)
    out.mkdir(parents=True, exist_ok=True)
    df.to_csv(out / "results.csv", index=False)
    summarize(df).to_csv(out / "summary.csv", index=False)
    (out / "summary.tex").write_text(to_latex(df, caption=caption))
    meta = {
        "qugrid_version": __version__,
        "numpy_version": np.__version__,
        "python": platform.python_version(),
        "platform": platform.platform(),
        "timestamp_utc": datetime.now(timezone.utc).isoformat(),
        "config": config or {},
    }
    (out / "config.json").write_text(json.dumps(meta, indent=2))
    return out

Classical references

Classical power flow: DC (linear) and Newton-Raphson AC.

Both follow the MATPOWER formulation, so results are directly comparable to what researchers get from rundcpf / runpf. Dense linear algebra keeps the code readable; the target scale (up to a few hundred buses) does not need sparsity.

DCResult dataclass

DC power flow solution.

Source code in src/qugrid/classical/power_flow.py

@dataclass
class DCResult:
    """DC power flow solution."""

    theta: np.ndarray  # bus voltage angles [rad]
    flow_mw: np.ndarray  # branch active flows, from-end [MW]
    slack_p_mw: float  # active power picked up by the slack bus [MW]
    loading: np.ndarray  # |flow| / RATE_A where a rating is given, else nan

    def max_loading(self) -> float:
        vals = self.loading[np.isfinite(self.loading)]
        return float(vals.max()) if vals.size else float("nan")

ACResult dataclass

Newton-Raphson AC power flow solution.

Source code in src/qugrid/classical/power_flow.py

@dataclass
class ACResult:
    """Newton-Raphson AC power flow solution."""

    vm: np.ndarray  # voltage magnitudes [pu]
    va: np.ndarray  # voltage angles [rad]
    converged: bool
    iterations: int
    mismatch_history: list[float] = field(default_factory=list)
    s_from_mva: np.ndarray | None = None  # complex from-end branch flows [MVA]
    s_to_mva: np.ndarray | None = None

    @property
    def v(self) -> np.ndarray:
        return self.vm * np.exp(1j * self.va)

    def losses_mw(self) -> float:
        return float(np.real(self.s_from_mva + self.s_to_mva).sum())

solve_dc(net)

Solve the DC power flow Bbus @ theta = P exactly.

This is also the linear system handed to quantum linear solvers (:func:qugrid.problems.dc_power_flow), so classical and quantum answers are comparable angle by angle.

Source code in src/qugrid/classical/power_flow.py

def solve_dc(net: Network) -> DCResult:
    """Solve the DC power flow ``Bbus @ theta = P`` exactly.

    This is also the linear system handed to quantum linear solvers
    (:func:`qugrid.problems.dc_power_flow`), so classical and quantum answers
    are comparable angle by angle.
    """
    bbus, bf, pbusinj, pfinj = net.bdc()
    pbus = np.real(net.sbus()) - net.bus[:, idx.GS] / net.baseMVA
    ref = net.ref
    keep = np.setdiff1d(np.arange(net.n_bus), [ref])
    theta = np.zeros(net.n_bus)
    theta[ref] = np.deg2rad(net.bus[ref, idx.VA])
    rhs = pbus[keep] - pbusinj[keep] - bbus[keep, ref] * theta[ref]
    theta[keep] = np.linalg.solve(bbus[np.ix_(keep, keep)], rhs)
    flow_pu = bf @ theta + pfinj
    flow_mw = flow_pu * net.baseMVA
    rate = net.branch[:, idx.RATE_A]
    loading = np.where(rate > 0, np.abs(flow_mw) / np.where(rate > 0, rate, 1.0), np.nan)
    slack_p = (bbus[ref] @ theta + pbusinj[ref]) * net.baseMVA + net.bus[ref, idx.PD]
    return DCResult(theta=theta, flow_mw=flow_mw, slack_p_mw=float(slack_p), loading=loading)

newton_raphson(net, tol=1e-08, max_iter=20)

Full Newton-Raphson AC power flow in polar coordinates.

Reference: MATPOWER's newtonpf (Zimmerman, Murillo-Sanchez, Thomas, "MATPOWER: Steady-State Operations, Planning and Analysis Tools for Power Systems Research and Education", IEEE Trans. Power Systems, 2011). Generator reactive-power limits are not enforced.

Source code in src/qugrid/classical/power_flow.py

def newton_raphson(net: Network, tol: float = 1e-8, max_iter: int = 20) -> ACResult:
    """Full Newton-Raphson AC power flow in polar coordinates.

    Reference: MATPOWER's ``newtonpf`` (Zimmerman, Murillo-Sanchez, Thomas,
    "MATPOWER: Steady-State Operations, Planning and Analysis Tools for Power
    Systems Research and Education", IEEE Trans. Power Systems, 2011).
    Generator reactive-power limits are not enforced.
    """
    ybus = net.ybus()
    sbus = net.sbus()
    pv, pq = net.pv, net.pq
    pvpq = np.concatenate([pv, pq]).astype(int)

    vm = net.bus[:, idx.VM].copy()
    va = np.deg2rad(net.bus[:, idx.VA])
    on = net.gen_on
    vm[net.gen_bus[on]] = net.gen[on, idx.VG]
    v = vm * np.exp(1j * va)

    history: list[float] = []
    converged = False
    it = 0
    for it in range(1, max_iter + 1):
        ibus = ybus @ v
        mis = v * np.conj(ibus) - sbus
        f_vec = np.concatenate([np.real(mis[pvpq]), np.imag(mis[pq])])
        norm = float(np.linalg.norm(f_vec, np.inf)) if f_vec.size else 0.0
        history.append(norm)
        if norm < tol:
            converged = True
            break

        diag_v = np.diag(v)
        diag_i = np.diag(ibus)
        diag_vn = np.diag(v / np.abs(v))
        ds_dva = 1j * diag_v @ np.conj(diag_i - ybus @ diag_v)
        ds_dvm = diag_v @ np.conj(ybus @ diag_vn) + np.conj(diag_i) @ diag_vn

        j11 = np.real(ds_dva[np.ix_(pvpq, pvpq)])
        j12 = np.real(ds_dvm[np.ix_(pvpq, pq)])
        j21 = np.imag(ds_dva[np.ix_(pq, pvpq)])
        j22 = np.imag(ds_dvm[np.ix_(pq, pq)])
        jac = np.block([[j11, j12], [j21, j22]])

        dx = np.linalg.solve(jac, -f_vec)
        n1 = len(pvpq)
        va[pvpq] += dx[:n1]
        vm[pq] += dx[n1:]
        v = vm * np.exp(1j * va)

    s_from, s_to = _branch_flows(net, v)
    return ACResult(
        vm=np.abs(v),
        va=np.angle(v),
        converged=converged,
        iterations=it,
        mismatch_history=history,
        s_from_mva=s_from,
        s_to_mva=s_to,
    )

Classical dispatch references: economic dispatch and enumerated unit commitment.

These provide the true optima against which QuGrid measures quantum solutions — including the error introduced by QUBO discretization itself, which is reported separately from solver error.

GenParams dataclass

One thermal generator for dispatch/commitment studies.

Cost is quadratic: c2 * P^2 + c1 * P + c0 [$ per hour], with c0 charged only when the unit is committed. startup is charged on each off-to-on transition.

Source code in src/qugrid/classical/dispatch.py

@dataclass
class GenParams:
    """One thermal generator for dispatch/commitment studies.

    Cost is quadratic: ``c2 * P^2 + c1 * P + c0`` [$ per hour], with ``c0``
    charged only when the unit is committed. ``startup`` is charged on each
    off-to-on transition.
    """

    name: str
    pmin: float
    pmax: float
    c2: float
    c1: float
    c0: float = 0.0
    startup: float = 0.0

    def cost(self, p: float | np.ndarray) -> float | np.ndarray:
        return self.c2 * p**2 + self.c1 * p + self.c0

UCSolution dataclass

A unit commitment schedule with its dispatch and cost.

Source code in src/qugrid/classical/dispatch.py

@dataclass
class UCSolution:
    """A unit commitment schedule with its dispatch and cost."""

    commit: np.ndarray  # (n_gen, n_periods) of 0/1
    power: np.ndarray  # (n_gen, n_periods) [MW]
    cost: float
    feasible: bool

economic_dispatch(gens, demand, tol=1e-09)

Exact single-period economic dispatch via bisection on marginal cost.

Solves min sum_i c_i(P_i) s.t. sum_i P_i = demand, pmin_i <= P_i <= pmax_i for strictly convex costs (c2 > 0). Returns (P, total_cost). Raises ValueError when the demand is outside the feasible range.

Source code in src/qugrid/classical/dispatch.py

def economic_dispatch(
    gens: list[GenParams], demand: float, tol: float = 1e-9
) -> tuple[np.ndarray, float]:
    """Exact single-period economic dispatch via bisection on marginal cost.

    Solves ``min sum_i c_i(P_i)  s.t.  sum_i P_i = demand, pmin_i <= P_i <= pmax_i``
    for strictly convex costs (``c2 > 0``). Returns ``(P, total_cost)``.
    Raises ``ValueError`` when the demand is outside the feasible range.
    """
    pmin = np.array([g.pmin for g in gens])
    pmax = np.array([g.pmax for g in gens])
    c1 = np.array([g.c1 for g in gens])
    c2 = np.array([g.c2 for g in gens])
    if np.any(c2 <= 0):
        raise ValueError("economic_dispatch requires strictly convex costs (c2 > 0)")
    if not (pmin.sum() - tol <= demand <= pmax.sum() + tol):
        raise ValueError(
            f"demand {demand:.6g} MW outside feasible range "
            f"[{pmin.sum():.6g}, {pmax.sum():.6g}] MW"
        )

    def output(lam: float) -> np.ndarray:
        return np.clip((lam - c1) / (2 * c2), pmin, pmax)

    lo = float(np.min(c1 + 2 * c2 * pmin))
    hi = float(np.max(c1 + 2 * c2 * pmax))
    for _ in range(200):
        lam = 0.5 * (lo + hi)
        total = output(lam).sum()
        if abs(total - demand) < tol:
            break
        if total < demand:
            lo = lam
        else:
            hi = lam
    p = output(0.5 * (lo + hi))
    # Absorb the residual rounding into an interior (unclamped) unit.
    residual = demand - p.sum()
    interior = np.flatnonzero((p > pmin + tol) & (p < pmax - tol))
    if interior.size and abs(residual) > 0:
        p[interior[0]] += residual
    cost = float(sum(g.cost(p[i]) for i, g in enumerate(gens)))
    return p, cost

solve_uc_enumerate(gens, demand, initial_on=None)

Exact unit commitment by enumerating all commitment patterns.

For each of the 2^(G*T) on/off patterns, the committed units are dispatched with exact economic dispatch. Intended for the small instances (G*T <= ~20) used to validate quantum solvers; it is the ground truth, not a production UC.

Source code in src/qugrid/classical/dispatch.py

def solve_uc_enumerate(
    gens: list[GenParams],
    demand: list[float] | np.ndarray,
    initial_on: np.ndarray | None = None,
) -> UCSolution:
    """Exact unit commitment by enumerating all commitment patterns.

    For each of the ``2^(G*T)`` on/off patterns, the committed units are
    dispatched with exact economic dispatch. Intended for the small instances
    (``G*T <= ~20``) used to validate quantum solvers; it is the ground truth,
    not a production UC.
    """
    demand = np.asarray(demand, dtype=float)
    n_gen, n_t = len(gens), len(demand)
    if n_gen * n_t > 24:
        raise ValueError("enumeration limited to G*T <= 24 variables")
    prev0 = np.zeros(n_gen) if initial_on is None else np.asarray(initial_on, dtype=float)

    best_cost = np.inf
    best: UCSolution | None = None
    for bits in product((0, 1), repeat=n_gen * n_t):
        commit = np.array(bits, dtype=float).reshape(n_gen, n_t)
        cost = 0.0
        power = np.zeros((n_gen, n_t))
        feasible = True
        for t in range(n_t):
            on = np.flatnonzero(commit[:, t])
            committed = [gens[i] for i in on]
            if not committed:
                feasible = False
                break
            try:
                p, c = economic_dispatch(committed, float(demand[t]))
            except ValueError:
                feasible = False
                break
            power[on, t] = p
            cost += c
            prev = prev0 if t == 0 else commit[:, t - 1]
            starts = np.maximum(commit[:, t] - prev, 0.0)
            cost += float(sum(g.startup * s for g, s in zip(gens, starts)))
        if feasible and cost < best_cost:
            best_cost = cost
            best = UCSolution(commit=commit, power=power, cost=cost, feasible=True)
    if best is None:
        raise ValueError("no feasible commitment pattern for the given demand")
    return best

Adapters

Bridges into the quantum ecosystems researchers already use.

QuGrid's built-in simulator is for zero-setup research; when you want shots, noise models, annealing hardware, or autodiff, export the same problem:

• :func:to_qiskit_operator / :func:to_qiskit_qaoa — IBM Qiskit
• :func:to_bqm / :func:result_from_sampleset — D-Wave Ocean (dimod)
• :func:to_pennylane_hamiltonian — Xanadu PennyLane

The same bridges also run in the other direction: :mod:external solvers <qugrid.adapters.external_solvers> hand the problem to a vendor solver and return the standard :class:~qugrid.solvers.base.Result, so qg.solve(problem, solver="dwave-sa") reads like any built-in solver.

Each function imports its ecosystem lazily and raises a clear message naming the extra to install (pip install qugrid[qiskit] and friends).

result_from_sampleset(problem, sampleset)

Wrap an Ocean SampleSet into a standard QuGrid :class:Result.

Source code in src/qugrid/adapters/dimod_adapter.py

def result_from_sampleset(problem: CombinatorialProblem, sampleset) -> Result:
    """Wrap an Ocean ``SampleSet`` into a standard QuGrid :class:`Result`."""
    _require_dimod()
    best = sampleset.first
    x = np.array([best.sample[i] for i in range(problem.qubo.n)], dtype=int)
    res = Result(solver=f"dimod:{type(sampleset).__name__}", problem=problem)
    finish_combinatorial(res, problem, x)
    res.resources["num_reads"] = len(sampleset)
    info = getattr(sampleset, "info", None)
    if info:
        timing = info.get("timing")
        if timing:
            res.resources["timing"] = timing
    return attach_reference(res, problem)

to_bqm(problem)

dimod.BinaryQuadraticModel of the problem's QUBO.

Submit it to any Ocean sampler::

from dwave.samplers import SimulatedAnnealingSampler
sampleset = SimulatedAnnealingSampler().sample(to_bqm(prob), num_reads=200)
result = result_from_sampleset(prob, sampleset)

Source code in src/qugrid/adapters/dimod_adapter.py

def to_bqm(problem: CombinatorialProblem):
    """``dimod.BinaryQuadraticModel`` of the problem's QUBO.

    Submit it to any Ocean sampler::

        from dwave.samplers import SimulatedAnnealingSampler
        sampleset = SimulatedAnnealingSampler().sample(to_bqm(prob), num_reads=200)
        result = result_from_sampleset(prob, sampleset)
    """
    _require_dimod()
    import dimod

    qubo = problem.qubo
    linear = {i: float(qubo.q[i, i]) for i in range(qubo.n) if qubo.q[i, i]}
    quadratic = {
        (i, j): float(2.0 * qubo.q[i, j])
        for i in range(qubo.n)
        for j in range(i + 1, qubo.n)
        if qubo.q[i, j]
    }
    return dimod.BinaryQuadraticModel(linear, quadratic, qubo.offset, vartype=dimod.BINARY)

solve_dimod_exact(problem, **_ignored)

Exhaustive enumeration with dimod.ExactSolver (n <= 24).

Source code in src/qugrid/adapters/external_solvers.py

def solve_dimod_exact(problem: CombinatorialProblem, **_ignored) -> Result:
    """Exhaustive enumeration with ``dimod.ExactSolver`` (n <= 24)."""
    _require("dimod", "dwave")
    import dimod

    from qugrid.adapters.dimod_adapter import to_bqm

    n = problem.qubo.n
    if n > MAX_EXACT_VARS:
        raise ValueError(f"refusing to enumerate 2^{n} states")
    res = Result(solver="dimod-exact", problem=problem)
    with timed(res.resources):
        sampleset = dimod.ExactSolver().sample(to_bqm(problem))
        x = _bits_from_sample(problem, sampleset.first.sample)
    finish_combinatorial(res, problem, x)
    res.resources.update({"n_vars": n, "states_enumerated": len(sampleset)})
    return attach_reference(res, problem)

solve_dwave_sa(problem, num_reads=100, num_sweeps=1000, seed=0, **_ignored)

Simulated annealing with dwave.samplers.SimulatedAnnealingSampler.

The Ocean stack's production annealer (C++ inner loop), run on the same QUBO as the built-in "sa" so the two are directly comparable.

Source code in src/qugrid/adapters/external_solvers.py

def solve_dwave_sa(
    problem: CombinatorialProblem,
    num_reads: int = 100,
    num_sweeps: int = 1000,
    seed: int = 0,
    **_ignored,
) -> Result:
    """Simulated annealing with ``dwave.samplers.SimulatedAnnealingSampler``.

    The Ocean stack's production annealer (C++ inner loop), run on the same
    QUBO as the built-in ``"sa"`` so the two are directly comparable.
    """
    _require("dwave.samplers", "dwave")
    from dwave.samplers import SimulatedAnnealingSampler

    from qugrid.adapters.dimod_adapter import to_bqm

    res = Result(solver="dwave-sa", problem=problem)
    with timed(res.resources):
        sampleset = SimulatedAnnealingSampler().sample(
            to_bqm(problem), num_reads=num_reads, num_sweeps=num_sweeps, seed=seed
        )
        x = _bits_from_sample(problem, sampleset.first.sample)
    finish_combinatorial(res, problem, x)
    res.resources.update(
        {
            "n_vars": problem.qubo.n,
            "num_reads": num_reads,
            "num_sweeps": num_sweeps,
            "seed": seed,
        }
    )
    return attach_reference(res, problem)

solve_qiskit_qaoa(problem, reps=2, maxiter=200, seed=0, **_ignored)

QAOA through qiskit-optimization's MinimumEigenOptimizer.

Shot-based sampling on Qiskit's statevector primitive, COBYLA outer loop. Slower than the built-in "qaoa" by two orders of magnitude at the same depth — the point is that a second implementation of the same algorithm lands on the same objective.

Source code in src/qugrid/adapters/external_solvers.py

def solve_qiskit_qaoa(
    problem: CombinatorialProblem,
    reps: int = 2,
    maxiter: int = 200,
    seed: int = 0,
    **_ignored,
) -> Result:
    """QAOA through qiskit-optimization's ``MinimumEigenOptimizer``.

    Shot-based sampling on Qiskit's statevector primitive, COBYLA outer
    loop. Slower than the built-in ``"qaoa"`` by two orders of magnitude at
    the same depth — the point is that a second implementation of the same
    algorithm lands on the same objective.
    """
    _require("qiskit_optimization", "qiskit")
    _require("qiskit_algorithms", "qiskit")
    from qiskit.primitives import StatevectorSampler
    from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager
    from qiskit_algorithms import QAOA
    from qiskit_algorithms.optimizers import COBYLA
    from qiskit_algorithms.utils import algorithm_globals
    from qiskit_optimization.algorithms import MinimumEigenOptimizer

    algorithm_globals.random_seed = seed  # fixes QAOA's random initial angles
    evaluations: list[int] = []
    res = Result(solver="qiskit-qaoa", problem=problem)
    with timed(res.resources):
        # The QAOA ansatz carries a PauliEvolutionGate. Synthesizing it once
        # through a pass manager, rather than at every objective call inside
        # the sampler, is worth two orders of magnitude in wall time.
        passes = generate_preset_pass_manager(
            optimization_level=1, basis_gates=["h", "rx", "rz", "cx"]
        )
        qaoa = QAOA(
            StatevectorSampler(seed=seed),
            COBYLA(maxiter=maxiter),
            reps=reps,
            transpiler=passes,
            callback=lambda count, *_rest: evaluations.append(count),
        )
        out = MinimumEigenOptimizer(qaoa).solve(_quadratic_program(problem))
        x = np.asarray(out.x, dtype=int)
    finish_combinatorial(res, problem, x)
    res.resources.update(
        {
            "n_qubits": problem.qubo.n,
            "p": reps,
            "evaluations": evaluations[-1] if evaluations else 0,
            "seed": seed,
        }
    )
    return attach_reference(res, problem)

to_pennylane_hamiltonian(problem)

Ising cost Hamiltonian as a PennyLane operator (autodiff-ready).

Use it in any PennyLane QNode, for example with qml.qaoa.cost_layer / qml.qaoa.mixer_layer.

Source code in src/qugrid/adapters/pennylane_adapter.py

def to_pennylane_hamiltonian(problem: CombinatorialProblem):
    """Ising cost Hamiltonian as a PennyLane operator (autodiff-ready).

    Use it in any PennyLane QNode, for example with
    ``qml.qaoa.cost_layer`` / ``qml.qaoa.mixer_layer``.
    """
    try:
        import pennylane as qml
    except ImportError as err:  # pragma: no cover
        raise ImportError(
            "pennylane is not installed; install the extra with "
            "`pip install qugrid[pennylane]`"
        ) from err

    ising = problem.qubo.to_ising()
    coeffs: list[float] = []
    ops: list = []
    if ising.offset:
        coeffs.append(float(ising.offset))
        ops.append(qml.Identity(0))
    for i in range(ising.n):
        if ising.h[i]:
            coeffs.append(float(ising.h[i]))
            ops.append(qml.PauliZ(i))
    for i in range(ising.n):
        for j in range(i + 1, ising.n):
            if ising.j[i, j]:
                coeffs.append(float(ising.j[i, j]))
                ops.append(qml.PauliZ(i) @ qml.PauliZ(j))
    return qml.Hamiltonian(coeffs, ops)

to_qiskit_operator(problem)

Ising cost Hamiltonian as a SparsePauliOp (identity term included).

Source code in src/qugrid/adapters/qiskit_adapter.py

def to_qiskit_operator(problem: CombinatorialProblem):
    """Ising cost Hamiltonian as a ``SparsePauliOp`` (identity term included)."""
    _require_qiskit()
    from qiskit.quantum_info import SparsePauliOp

    ising = problem.qubo.to_ising()
    n = ising.n
    terms: list[tuple[str, list[int], float]] = []
    if ising.offset:
        terms.append(("I", [0], float(ising.offset)))
    for i in range(n):
        if ising.h[i]:
            terms.append(("Z", [i], float(ising.h[i])))
    for i in range(n):
        for j in range(i + 1, n):
            if ising.j[i, j]:
                terms.append(("ZZ", [i, j], float(ising.j[i, j])))
    return SparsePauliOp.from_sparse_list(terms, num_qubits=n)

to_qiskit_qaoa(problem, p=2)

Parameterized QAOA circuit for the problem.

Returns (circuit, gammas, betas) where the parameter vectors bind the p cost and mixer angles. Ready for Qiskit primitives, transpilers, and hardware backends.

Source code in src/qugrid/adapters/qiskit_adapter.py

def to_qiskit_qaoa(problem: CombinatorialProblem, p: int = 2):
    """Parameterized QAOA circuit for the problem.

    Returns ``(circuit, gammas, betas)`` where the parameter vectors bind the
    ``p`` cost and mixer angles. Ready for Qiskit primitives, transpilers,
    and hardware backends.
    """
    _require_qiskit()
    from qiskit import QuantumCircuit
    from qiskit.circuit import ParameterVector

    ising = problem.qubo.to_ising()
    n = ising.n
    gammas = ParameterVector("gamma", p)
    betas = ParameterVector("beta", p)
    qc = QuantumCircuit(n)
    qc.h(range(n))
    pairs = [(i, j) for i in range(n) for j in range(i + 1, n) if ising.j[i, j]]
    for layer in range(p):
        for i in range(n):
            if ising.h[i]:
                qc.rz(2.0 * float(ising.h[i]) * gammas[layer], i)
        for i, j in pairs:
            qc.rzz(2.0 * float(ising.j[i, j]) * gammas[layer], i, j)
        qc.rx(2.0 * betas[layer], range(n))
    qc.measure_all()
    return qc, gammas, betas