Diagnostics & Validation

How to detect the failure modes from Bounds, Branches, and Feasibility, rather than merely avoid them. This page is aimed at the validation audience: you are checking your own linear, convex, or nonconvex formulation against the benchmark and need to know what to measure and what a discrepancy means.

The throughline: two well-known generic optimization tools become your two domain-specific diagnostics. Bound everything and see what binds is the unboundedness test of §1; evaluating a relaxed solution against the AC equations is the exactness test of §5.

BMOPFTools-native shortcuts

The JuMP recipes below apply to your model. Two of these diagnostics already exist inside this package, and you should reach for them first:

AC-feasibility of a candidate point (Symptom 3) can be tested by first pinning its decisions with project_solution, then running solve_feasibility_opf on that determined snapshot. A converged near-zero slack plus independently checked residuals supplies a numerical feasible point. This is the native, four-wire analogue of primal_feasibility_report against the AC model.
Ill-posedness before you solve — flat objectives, free injections, and cost degeneracy — is caught structurally by the benchmark-readiness flags (W.BENCH.GEN_ZERO_COST, W.BENCH.GEN_DEGENERATE_COST, W.BENCH.GEN_NO_DOF; see Methodology § Benchmark readiness). Diverging iterates (NORM_LIMIT) from a zero-cost generator are the runtime echo of W.BENCH.GEN_ZERO_COST.

Symptom 1 — `INFEASIBLE_OR_UNBOUNDED` / `DUAL_INFEASIBLE`

First, disambiguate. Unboundedness can only come from the objective, so drop it and re-solve (YALMIP; JuMP):

@objective(model, Min, 0)   # feasibility test
optimize!(model)
# feasible now  ⇒ original was UNBOUNDED
# infeasible now ⇒ original was (also) INFEASIBLE

Then localize the unboundedness. Add large finite bounds to every free or one-sided variable, re-solve, and inspect which variables sit at the artificial bound — those are your unbounded directions:

for v in all_variables(model)
    has_lower_bound(v) || set_lower_bound(v, -1e6)
    has_upper_bound(v) || set_upper_bound(v,  1e6)
end
optimize!(model)
for v in all_variables(model)
    if isapprox(value(v), 1e6; atol=1) || isapprox(value(v), -1e6; atol=1)
        @info "at artificial bound" v value(v)
    end
end

In a power-flow model the variables pinned at the artificial bound are almost always voltage magnitudes (missing $\underline v$ / $\overline v$) or an unbounded generator — i.e. exactly the missing operational limits from §1. Erwin Kalvelagen's variant — bind the objective to a large-bounded auxiliary variable via an equality, then read its value — is the same idea (the canonical infeasible/unbounded writeup).

Symptom 2 — model is infeasible and you do not know why

Use an irreducible infeasible subsystem (IIS) if your solver supports it:

compute_conflict!(model)
if get_attribute(model, MOI.ConflictStatus()) == MOI.CONFLICT_FOUND
    iis, _ = copy_conflict(model)
    print(iis)
end

When no IIS is available, the penalty relaxation locates the binding constraints by letting them be violated at a cost:

map = relax_with_penalty!(model)   # adds slacks + penalty to the objective
optimize!(model)
for (con, slack) in map
    value(slack) > 1e-6 && @info "violated" con value(slack)
end

Distinguish modeling infeasibility from physical infeasibility

relax_with_penalty! does not relax variable bounds or integrality. More important here: in a constant-power model an infeasibility may be physical — you are past the loadability / collapse boundary (§4), not mis-modeled. If the IIS centres on the power-balance and voltage-bound constraints of a heavily loaded sub-network, suspect collapse, and confirm with a continuation sweep rather than editing constraints.

Symptom 3 — is my relaxation exact? (the central validation check)

Evaluate the relaxed solution against the nonconvex AC constraints and read off the violation. JuMP's primal_feasibility_report does this against any model:

# Schematic — solve the relaxation, then test its point against the AC model.
# NOTE: the two models do not share a variable space (ℓ, w vs V, S), so you must
# first map the relaxed solution into the AC model's variables (recover V from w and
# the flows, keyed by name/bus), not copy variable objects across models:
point  = recover_ac_point(relaxed_model)          # your relaxation → AC-variable map
report = primal_feasibility_report(ac_model, point)
isempty(report) ? @info("AC-feasible: relaxation exact") :
                  @info("AC-infeasible: relaxation INEXACT", report)

In BMOPFTools the same check is available natively for a four-wire network: pin the candidate dispatch into the case and call solve_feasibility_opf, which minimises an injected slack current. A converged near-zero slack, together with independently recomputed residuals, demonstrates a valid numerical power flow. A non-zero slack localises the residual used by that relaxed solve:

snap    = project_solution(net, result)
res     = solve_feasibility_opf(snap; optimizer = Ipopt.Optimizer)
slack_A = res["total_slack_magnitude_A"]
slack_A < 1e-3 ? @info("candidate has near-zero KCL slack; verify residuals") :
                 @info("candidate/relaxation has residual current", slack_A)

For the branch-flow relaxation on a radial network you can also read the gap directly per branch:

\[\text{gap}_{ij} \;=\; \ell_{ij} \;-\; \frac{|S_{ij}|^2}{v_i}\,,\]

which is $\approx 0$ on every branch when exact and strictly positive where the cone is slack. For the SDP / BIM relaxation, check the per-clique voltage matrix rank (or the ratio of the two largest eigenvalues): a rank-1 matrix (eigenratio $\to \infty$, second eigenvalue $\to 0$) means exact (Lupien & Lesage-Landry, 2023).

Cone tightness is not sufficient on meshed networks

On a mesh, all branch cones can bind while the point is still AC-infeasible: the SOC branch-flow relaxation drops the cycle angle-consistency (KVL-around-loops) constraints, so a solution with every $\text{gap}_{ij}\approx 0$ can fail to admit any consistent voltage-angle assignment around a loop. Per-branch tightness certifies exactness only on radial feeders. On meshes, fall back to the direct AC-residual checks above — primal_feasibility_report against the AC model, or solve_feasibility_opf's slack current — which remain valid regardless of topology because they test the recovered point against the full AC equations rather than the relaxation's own constraints.

Where to expect inexactness

Empirically the relaxation goes inexact precisely when upper voltage bounds bind and their duals exceed a threshold — and exact solutions reliably show binding upper bounds on active/reactive withdrawals instead (Gan et al., 2015; Bobo et al., 2020). If your validation diverges from the benchmark, check whether you are in a loss-rewarding objective (list) or against an upper voltage bound first.

Symptom 4 — I think I am on the wrong (low-voltage) branch

A nonconvex solver can converge to the low-voltage solution, especially under a loss-rewarding objective or a poor start. Checks:

Voltage profile. Operational solutions often sit near $1$ p.u.; a solution with a cluster of buses well below $\underline v$-class values (e.g. $0.5$–$0.8$ p.u.) on a normally loaded feeder is a warning signal. Magnitude alone does not classify a voltage sheet. Confirm the branch under a declared load/control continuation using the power-flow Jacobian and voltage response.
Restart from a high-voltage start. Flat start or warm-start from a linear-model solution. Under the assumptions of Dvijotham et al. the standard fixed-point methods recover the high-voltage solution (2017); do not export that theorem to an arbitrary meshed, multiphase OPF. Distinct verified terminal states establish multiplicity. Agreement across tested starts is evidence, not proof of uniqueness.
Cross-solver / cross-formulation. Solve an equivalent formulation or relaxation too. If a relaxation is exact, it supplies an AC-feasible global optimum, not necessarily a unique or upper-sheet optimum. Check the recovered state independently and compare only models with identical objective, bounds, angle conventions, and branch-rating semantics (§2).

Symptom 5 — the solve is slow, churns in restoration, or the duals look wrong

A correct model can still misbehave numerically. The triage is to decide which of three causes you are looking at — physical, degeneracy, or non-smoothness — before touching the formulation (Trusting the solver §4–§6):

Physical (collapse). A heavily loaded sub-network past the nose (§4). The IIS / penalty relaxation from Symptom 2 centres on the power-balance and voltage-bound constraints there. Confirm with a continuation sweep, not by editing constraints.
Degeneracy (CQ failure). Slow convergence with a stable primal but unstable duals is the signature (Trusting the solver §4). Look for pinned or linearly dependent active constraints — most often the W.BENCH.GEN_* structures (Methodology). Independently verify the primal state; treat local-optimality, price, and sensitivity claims as unsupported until the CQ failure is resolved.
Non-smoothness / zero voltage. Restoration churn with exploding Hessian entries and vanishing steps points at a fractional ZIP exponent or a zero-voltage bilinear bifurcation (Trusting the solver §5–§6). Check that a strictly positive voltage floor and the load-model floor are in place.

Two cheap cross-checks separate "the solver struggled" from "the answer is wrong":

# 1. Multi-start: does the verdict / objective survive a change of start?
starts = [flat_start, highvoltage_start, warmstart_from_linear]
for s in starts
    set_start_value.(all_variables(model), s)
    optimize!(model)
    @info "start" termination_status(model) objective_value(model)
end
# divergent verdicts ⇒ start-dependence (multiplicity or numerical), not a fixed answer

# 2. Conditioning of the active-constraint Jacobian (LICQ proxy):
#    a tiny smallest singular value flags a constraint-qualification failure.
J = jacobian_of_active_constraints(model)   # your assembly of ∇g for active g
σ = svdvals(Matrix(J))
@info "active-Jacobian conditioning" σ_min=minimum(σ) σ_max=maximum(σ)
# σ_min ≈ 0 ⇒ degenerate active set (Trusting the solver §4)

For the AC-feasibility leg of this triage, reach for solve_feasibility_opf: a verified near-zero total_slack_magnitude_A supplies a feasible point, while non-zero slack is a local diagnostic that should be checked across starts and formulations.

Validation checklist against the benchmarks

Reproduce the benchmark's reported objective with your formulation and the same objective and bounds.
If you use a relaxation, run Symptom 3 and confirm exactness; only compare the value if exact.
If values differ, classify before debugging code: is the row in the decision matrix a ✗ cell? Is the objective on the loss-maximization list? Are upper voltage bounds binding? Expected divergence is not a bug.
For collapse/loadability instances, compare margins from a continuation sweep, not single-point objective values.