Checking the DAE index#

drto.check_index tests whether a declared model is an index-one DAE: whether every algebraic variable is determined, pointwise, by the algebraic constraints. A model that fails is higher-index, and its symptoms in optimization are indirect — undetermined variables at the first time point, nearly dependent constraint rows, stalled dual iterations — so the check names the offending variables and constraints up front instead.

The check runs in two layers. The structural layer matches the algebraic constraints to the algebraic variables on the incidence graph: a variable no constraint can pair with makes the model structurally higher-index, and Pantelides’ algorithm then reports the structural index. The numerical layer evaluates the algebra’s Jacobian at the model’s current values, equilibrates it so units drop out, and estimates its condition, since a structurally sound algebra can still be singular at the point.

The pendulum, four ways#

The classic index ladder (APMonitor’s pendulum): a mass on a cord of length \(s\), positions \(x, y\), velocities \(v, w\), and cord tension \(\lambda\). The equations of motion are the same four differential equations in every version,

\[\dot x = v, \quad \dot y = w, \quad m\dot v = -2x\lambda, \quad m\dot w = -mg - 2y\lambda,\]

and the versions differ in the fifth equation, the one that must determine \(\lambda\):

index 3 — the cord length: \(x^2 + y^2 = s^2\). No \(\lambda\) anywhere in the algebra; two differentiations stand between the constraint and the tension.
index 2 — the length constraint differentiated once: \(xv + yw = 0\). Still no \(\lambda\).
index 1 — differentiated twice: \(m(v^2 + w^2 - gy) - 2\lambda(x^2 + y^2) = 0\). The tension appears, and the algebra determines it.
index 0 — differentiated a third time, giving \(\lambda\) its own differential equation: the model is an ODE and \(\lambda\) is a fifth state.

import pyomo.environ as pyo
from pyomo.dae import ContinuousSet, DerivativeVar

import drto

G, MASS, CORD = 9.81, 1.0, 1.0


def pendulum(index):
    m = pyo.ConcreteModel()
    m.time = ContinuousSet(bounds=(0, 1))
    m.x = pyo.Var(m.time, initialize=0.5)
    m.y = pyo.Var(m.time, initialize=-0.866)
    m.v = pyo.Var(m.time, initialize=0.1)
    m.w = pyo.Var(m.time, initialize=0.1)
    m.lam = pyo.Var(m.time, initialize=0.1)
    m.dx = DerivativeVar(m.x, wrt=m.time, initialize=0)
    m.dy = DerivativeVar(m.y, wrt=m.time, initialize=0)
    m.dv = DerivativeVar(m.v, wrt=m.time, initialize=0)
    m.dw = DerivativeVar(m.w, wrt=m.time, initialize=0)
    if index == 0:
        m.dlam = DerivativeVar(m.lam, wrt=m.time, initialize=0)

    drto.horizon(m.time)
    pyo.TransformationFactory("dae.collocation").apply_to(
        m, wrt=m.time, nfe=2, ncp=2, scheme="LAGRANGE-RADAU")

    m.ode_x = pyo.Constraint(m.time, rule=lambda b, t: b.dx[t] == b.v[t])
    m.ode_y = pyo.Constraint(m.time, rule=lambda b, t: b.dy[t] == b.w[t])
    m.ode_v = pyo.Constraint(
        m.time, rule=lambda b, t: MASS * b.dv[t] == -2 * b.x[t] * b.lam[t])
    m.ode_w = pyo.Constraint(
        m.time, rule=lambda b, t:
        MASS * b.dw[t] == -MASS * G - 2 * b.y[t] * b.lam[t])
    dyn = [m.ode_x, m.ode_y, m.ode_v, m.ode_w]
    if index == 3:
        m.alg = pyo.Constraint(
            m.time, rule=lambda b, t: b.x[t] ** 2 + b.y[t] ** 2 == CORD ** 2)
    elif index == 2:
        m.alg = pyo.Constraint(
            m.time, rule=lambda b, t:
            b.x[t] * b.v[t] + b.y[t] * b.w[t] == 0)
    elif index == 1:
        m.alg = pyo.Constraint(
            m.time, rule=lambda b, t:
            MASS * (b.v[t] ** 2 + b.w[t] ** 2 - G * b.y[t])
            - 2 * b.lam[t] * (b.x[t] ** 2 + b.y[t] ** 2) == 0)
    else:
        m.ode_lam = pyo.Constraint(
            m.time, rule=lambda b, t:
            b.dlam[t] * (b.x[t] ** 2 + b.y[t] ** 2)
            == -4 * b.lam[t] * (b.x[t] * b.v[t] + b.y[t] * b.w[t])
            - 1.5 * G * MASS * b.w[t])
        dyn.append(m.ode_lam)

    states = [m.x, m.y, m.v, m.w] + ([m.lam] if index == 0 else [])
    drto.state(*states)
    drto.dynamics(*dyn)
    return m

The check reads each version the way the ladder predicts: the index-3 and index-2 forms fail structurally with the tension named as the variable nothing determines, and Pantelides reports the depth; the index-1 form passes with a clean condition estimate; the ODE form passes with the algebra empty.

for k in (3, 2, 1, 0):
    print(f"pendulum written as index {k}:")
    print(drto.check_index(pendulum(k)))
    print()

pendulum written as index 3:
drto check_index at t = 0.166667
  algebra: 1 constraints, 1 algebraic variables
  structural: FAILED - variables no algebraic constraint determines:
    lam[0.166667]
  unmatched constraints:
    alg[0.166667]
  structural index: 3, computed from which variables appear in which equations (Pantelides' algorithm); a coefficient that cancels numerically can hide a dependency the pattern shows, so the true index can be higher
  note: Dulmage-Mendelsohn: 1 constraints in the overconstrained subsystem, 1 variables in the underconstrained subsystem
  verdict: not index one: structurally higher index, see the unmatched variables

pendulum written as index 2:
drto check_index at t = 0.166667
  algebra: 1 constraints, 1 algebraic variables
  structural: FAILED - variables no algebraic constraint determines:
    lam[0.166667]
  unmatched constraints:
    alg[0.166667]
  structural index: 2, computed from which variables appear in which equations (Pantelides' algorithm); a coefficient that cancels numerically can hide a dependency the pattern shows, so the true index can be higher
  note: Dulmage-Mendelsohn: 1 constraints in the overconstrained subsystem, 1 variables in the underconstrained subsystem
  verdict: not index one: structurally higher index, see the unmatched variables

pendulum written as index 1:
drto check_index at t = 0.166667
  algebra: 1 constraints, 1 algebraic variables
  structural: full matching
  numerical: condition estimate 1.000e+00
  verdict: index one

pendulum written as index 0:
drto check_index at t = 0.166667
  algebra: 0 constraints, 0 algebraic variables
  structural: full matching
  numerical: nothing to evaluate
  verdict: index zero: the pointwise algebra is empty, the model is an ODE

The solvent extraction stage, two ways#

The same physics, posed twice — the case that motivated the check. The MSContactor form (models/prommis_sx.py) keeps PrOMMiS’s transfer extents: variables that appear only in the balances, with the distribution equilibria written on the concentrations. No algebraic constraint determines an extent — exactly the pendulum’s tension, a hundred times over — and every symptom of a higher-index model followed it through this example’s development: extents undetermined at the first time point, nearly dependent rows, dual iterations that stall at every scaling.

The check says all of that in one report. The pointwise algebra is not even square, the extents head the unmatched variables, and the distribution constraints head the unmatched rows.

from models import prommis_sx

m_extent = prommis_sx.build(N=1, h=0.25)
print(drto.check_index(m_extent))

drto check_index at t = 0.25
  algebra: 514 constraints, 439 algebraic variables
  structural: FAILED - variables no algebraic constraint determines:
    fs.ms.aqueous_settler[1].unit.properties[0.25,1.0].flow_vol
    fs.ms.aqueous_settler[1].unit.properties[0.25,1.0].conc_mass_comp[H2O]
    fs.ms.aqueous_settler[1].unit.properties[0.25,1.0].conc_mass_comp[HSO4]
    fs.ms.organic_settler[1].unit.properties[0.25,1.0].flow_vol
    fs.ms.organic_settler[1].unit.properties[0.25,1.0].conc_mass_comp[Kerosene]
    fs.ms.mixer[1].unit.mscontactor.heterogeneous_reaction_extent[0.25,1,La_mass_transfer]
    fs.ms.mixer[1].unit.mscontactor.heterogeneous_reaction_extent[0.25,1,Y_mass_transfer]
    fs.ms.mixer[1].unit.mscontactor.heterogeneous_reaction_extent[0.25,1,Pr_mass_transfer]
    fs.ms.mixer[1].unit.mscontactor.heterogeneous_reaction_extent[0.25,1,Ce_mass_transfer]
    fs.ms.mixer[1].unit.mscontactor.heterogeneous_reaction_extent[0.25,1,Nd_mass_transfer]
    ... and 30 more
  unmatched constraints:
    fs.ms.mixer[1].unit.distribution_extent_constraint[0.25,1,La]
    fs.ms.mixer[1].unit.distribution_extent_constraint[0.25,1,Y]
    fs.ms.mixer[1].unit.distribution_extent_constraint[0.25,1,Pr]
    fs.ms.mixer[1].unit.distribution_extent_constraint[0.25,1,Ce]
    fs.ms.mixer[1].unit.distribution_extent_constraint[0.25,1,Nd]
    fs.ms.mixer[1].unit.distribution_extent_constraint[0.25,1,Sm]
    fs.ms.mixer[1].unit.distribution_extent_constraint[0.25,1,Gd]
    fs.ms.mixer[1].unit.distribution_extent_constraint[0.25,1,Dy]
    fs.ms.mixer[1].unit.distribution_extent_constraint[0.25,1,Al]
    fs.ms.mixer[1].unit.distribution_extent_constraint[0.25,1,Ca]
    ... and 105 more
  structural index: at least 2 (the pointwise algebra cannot determine its variables, and differentiation cannot balance an overdetermined pointwise system)
  note: 635 constraints without the time coordinate are outside the pointwise algebra
  note: 918 variables without the time coordinate are outside the pointwise algebra
  note: the pointwise algebra is not square: 514 constraints for 439 algebraic variables
  note: Dulmage-Mendelsohn: 117 constraints in the overconstrained subsystem, 113 variables in the underconstrained subsystem
  verdict: not index one: structurally higher index, see the unmatched variables

The reformulated stage (models/prommis_sx_index_one.py) poses the balances on the reaction invariants — each metal’s total across both phases, the hydrogen, sulfur, and extractant combinations — so the transfer extents never exist and the equilibria split each total algebraically. The check confirms the posing: full matching, and a condition estimate at order 1e3 after equilibration.

from models import prommis_sx_index_one

m_invariant = prommis_sx_index_one.build(N=1, h=0.25)
print(drto.check_index(m_invariant))

drto check_index at t = 0.25
  algebra: 215 constraints, 215 algebraic variables
  structural: full matching
  numerical: condition estimate 6.510e+03
  note: 131 constraints without the time coordinate are outside the pointwise algebra
  note: 2 variables without the time coordinate are outside the pointwise algebra
  verdict: index one

One report, before anything solves: the extent form is structurally higher-index with the culprits named, and the invariant form is index one. The three weeks of symptoms between those two sentences are what the check exists to skip.

Checking the DAE index

Contents

Checking the DAE index#

The pendulum, four ways#

The solvent extraction stage, two ways#