Impasse Point Algorithm Implementation

Implementation of the Rabier & Rheinboldt (1994) algorithm for integrating differential-algebraic equations (DAEs) through impasse points (singularities).

Reference

P.J. Rabier and W.C. Rheinboldt, "On the Computation of Impasse Points of Quasi-Linear Differential-Algebraic Equations," Mathematics of Computation, Vol. 62, No. 205, pp. 133-154, January 1994. DOI: 10.1090/S0025-5718-1994-1208224-6

File Structure

• impasse_algorithm.py: Core library implementing the algorithm from the paper
• run_Rabier_Rheinboldt_1994_circuit.py: Circuit example from the original paper (working)
• run_thesis_circuit.py: Diode circuit from thesis (has scaling issues)

Current Implementation Status

What Works ✅

1. Core Algorithm: Full implementation of the Rabier & Rheinboldt (1994) algorithm including:

   • Constraint manifold projection using Newton iteration
   • Normal and tangent space basis computation via QR factorization
   • Augmented system formulation for singularity traversal
   • Adaptive Runge-Kutta integration (DOPRI5) with error control
   • Impasse point detection via gamma sign changes

2. Paper Circuit Example: The circuit from Rabier & Rheinboldt (1994) successfully reproduces the paper's results:

   • Integrates through impasse points
   • Handles orientation flip after crossing singularity
   • Matches published trajectories

3. IDA Comparison: Integration using SUNDIALS IDA solver for validation

Known Issues ⚠️

1. Thesis Diode Circuit: The autonomous diode circuit has severe scaling problems:

   • Matrix Bc = A1 @ Uc has extremely small entries (~1e-11)
   • Gamma becomes nearly zero (~1e-21), causing integration to stall
   • Real time t doesn't advance (stuck at ~1e-15 while pseudo-time s progresses)
   • Root cause: Different scales in A1 matrix (inductor equation has L=0.01, time equation has 1)
   • Note: This system does not necessarily have an impasse point, but more than once in the paper the authors claimed the algorithm could be used to solve through other kinds of singularities

2. Exponential Overflow Protection: Diode model uses piecewise linearization:

   • For v_D/V_th < 60: exponential model
   • For v_D/V_th ≥ 60: linear continuation
   • Uses np.where for both scalar and array inputs

Open Questions in Code Comments

Several questions remain about the algorithm's behavior (see code comments):

• Proper handling of gamma sign changes
• Orientation flipping strategies
• Relationship between pseudo-time and real time at impasse points

DAE System Definition

The code uses an abstract base class for defining DAE systems:

class DAESystem(ABC):
    """
    DAE system of the form:
        A1(x) @ x' = G1(x)    (r differential equations)
        G2(x) = 0             (p algebraic constraints)

    Dimensions:
        n = total variables
        p = number of constraints
        r = n - p (differential equations)
    """

    @abstractmethod
    def A1(self, x):
        """Return differential equation matrix (r x n)"""

    @abstractmethod
    def G1(self, x):
        """Return differential equation RHS (r,)"""

    @abstractmethod
    def G2(self, x):
        """Return algebraic constraints (p,)"""

    @abstractmethod
    def DG2(self, x):
        """Return Jacobian of G2 (p x n)"""

Usage

Basic Integration

import impasse_algorithm as imp

# Create system
system = imp.CircuitRabierRheinboldt1994(gamma_param=1.0)

# Get initial condition
x0 = imp.get_initial_condition("paper", system=system)

# Integrate through impasse point
sol, time = imp.integrate(
    x0,
    t_final=1.5,
    ds0=0.0001,              # Initial pseudo-time step
    system=system,
    flip_after_crossing=True, # Flip orientation at singularity
    initial_flip=True,        # Start with flipped orientation
    rk_method=imp.DOPRI5,
    verbose=False
)

Running Examples

Paper circuit (working):

python run_Rabier_Rheinboldt_1994_circuit.py

Thesis circuit (has issues):

python run_thesis_circuit.py

Integration Methods

Impasse point algorithm (for systems with singularities):

sol, time = imp.integrate(
    x0,
    t_final=1.5,
    ds0=0.0001,
    system=system,
    max_ds=1.0,               # Maximum pseudo-time step
    reltol=1e-6,              # Relative tolerance
    abstol=1e-8,              # Absolute tolerance
    flip_after_crossing=True,
    initial_flip=False,
    rk_method=imp.DOPRI5,
    verbose=True              # Print progress
)

IDA solver (standard DAE solver, for comparison): For the IDA solver, installing SUNDIALS is required. It is available on Homebrew.

sol, time = imp.integrate_IDA(
    x0,
    t_final=1.5,
    system=system,
    xdot0=None,                    # Initial derivatives (computed if None)
    dt=1e-6,                       # Time step
    algebraic_vars_idx=[0, 1, 2]  # Indices of algebraic variables
)

Algorithm Details

Key Features

1. Pseudo-time Parametrization: Uses pseudo-time s instead of real time t:

   • dt/ds = gamma(y(s)) where gamma can approach zero at impasse points
   • Allows integration to continue when standard solvers fail

2. Adaptive Step Size: DOPRI5 embedded Runge-Kutta with error control

3. Constraint Projection: Newton iteration projects solution onto constraint manifold

4. Singularity Detection: Monitors gamma sign changes to detect impasse point crossings

Verbose Output

Setting verbose=True provides integration progress:

Starting integration: t_final=1.5, ds0=0.0001
Initial state: x0=[...]
Step 100: t=0.082847, ds=0.000987, gamma=0.973421, crossings=0
Singularity crossing detected at step 815, t=0.798234
Orientation flipped
Integration complete: 1523 steps, 1 crossings
Final time: 1.500123 (target: 1.5)

Limitations and Future Work

1. Scaling Sensitivity: Algorithm can fail when system matrices have vastly different scales (see thesis circuit)

2. Physical Interpretation: Gamma sign changes and orientation flipping need careful interpretation for physical systems

3. Autonomous Systems: Time must be added as a variable; proper handling of time scaling is important

Dependencies

• numpy: Array operations
• scipy: Linear algebra (SVD, QR factorization)
• matplotlib: Plotting
• scikits.odes: IDA solver (SUNDIALS wrapper... so installing SUNDIALS is required)

Creating Custom Systems

To implement a new DAE system:

1. Inherit from DAESystem
2. Implement A1, G1, G2, DG2 methods
3. Ensure initial condition satisfies constraints

See CircuitRabierRheinboldt1994 for a complete example.