Class DaeAnalyzer

All Implemented Interfaces:
FluidAnalyzer

public class DaeAnalyzer extends ClosingAndStateDepMethodsAnalyzer
Differential-algebraic formulation of the min-normal closure, backing options.method='dae'. Java twin of MATLAB's solver_fluid_dae.

MinNormalAnalyzer already solves a differential system (the mean) coupled to an algebraic one (the covariance). It solves them by SUCCESSIVE SUBSTITUTION: integrate the mean to its fixed point at a held variance, solve the Lyapunov equation there, extract sigma2, repeat, up to 20 times and only to CoarseTol. This analyzer states the same closure as one system and solves it as one system. NOTHING ABOUT THE CLOSURE CHANGES -- the drift, the rate factors and the Lyapunov equation are taken unmodified from FluidMomentTerms and FluidLyapunov -- only the way the coupled equations are discharged.

Two modes, chosen by the horizon:

STEADY STATE (unbounded horizon, the usual case) solves

     0 = D r(x, sigma2)               drift residual, one row per state
     0 = C x - Nchain                 population conservation, per closed chain
     0 = sigma2 - sigmaOf(x, sigma2)  closure consistency, per closable station
 

simultaneously by a damped projected Newton. THE FIXED POINT IS NOT FOUND BY INTEGRATING TO IT: integrating a stable ODE until it stops moving is a poor way to solve f(x)=0, because the cost is set by the slowest mode of the model rather than by the accuracy wanted. One seed trajectory is still integrated, cheaply and at the first-order closure, because Newton needs a point inside the basin; everything after that is algebraic.

TRANSIENT (finite horizon) integrates the index-1 DAE with a SINGULAR mass matrix, the covariance advancing alongside the mean. This and kp are the only fluid methods producing a time-varying second moment; minnormal evaluates its whole transient at the single stationary variance. The integrator is Rodas, because LSODA solves y' = f and cannot carry a mass matrix at all.

WHAT THE ALGEBRAIC CONSTRAINT BUYS. Population conservation otherwise holds only to integrator tolerance: it is a CONSEQUENCE of the drift (the rows of D sum to zero on a closed chain), never an equation. Writing it as a constraint enforces it to solver tolerance, and it is also what makes the Newton system solvable at all -- the drift Jacobian is singular along exactly the conserved directions, the same singularity FluidLyapunov works around by projecting onto range(D), so the constraint rows supply the missing rank instead of a pseudo-inverse hiding it.

WHY THE COVARIANCE IS NOT A NEWTON UNKNOWN. Sigma is nstate^2 entries, so a Jacobian over it is quartic work -- strictly worse than the cubic Lyapunov solves it would replace. Sigma is LINEAR in itself for a held x, so it is eliminated by one Lyapunov solve per residual evaluation and only sigma2, M numbers, joins x in the unknown vector.

Java 8: no var, no List.of, no switch expressions.

See Also:
  • Field Details

    • sigmaMatrix

      public Matrix sigmaMatrix
      State-level stationary covariance at the converged fixed point.
    • qVar

      public Matrix qVar
      Per-(station,class) queue-length variance.
    • classBlock

      public int[][] classBlock
      State coordinates of each (station,class), flattened as i*K+k.
    • stationBlock

      public int[][] stationBlock
      State coordinates of each station.
    • sigma2

      public double[] sigma2
      Per-station population variance.
    • sigma2Drift

      public double[] sigma2Drift
      The same variance as it enters the DRIFT: zero where nothing is closed.
    • outerIters

      public int outerIters
      Newton steps taken, summed over the active-set passes.
    • residual

      public double residual
      Infinity norm of the residual the solve stopped at.
    • converged

      public boolean converged
      Whether that residual reached options.tol.
    • conservation

      public double conservation
      Largest violation of the conservation rows at the reported point.
    • sigmat

      public Matrix[] sigmat
      Transient state covariance, one matrix per reported time.
    • qVart

      public Matrix[] qVart
      Transient per-(station,class) variance, flattened as i*K+k.
    • tvar

      public Matrix tvar
      The time grid the transient covariance is reported on.
    • capacityLabel

      public String[] capacityLabel
      Human-readable origin of each capacity constraint.
    • capacityB

      public Matrix capacityB
      Right-hand side of each capacity constraint.
    • capacityValue

      public Matrix capacityValue
      Value each constrained quantity actually took.
    • capacityActive

      public int[] capacityActive
      Which constraints bound at the reported point.
    • staging

      public Matrix staging
      Mass held OUTSIDE a capped region, per staging coordinate.
    • stagingRegion

      public int[] stagingRegion
      Region each staging coordinate belongs to.
    • stagingClass

      public int[] stagingClass
      Class each staging coordinate carries.
    • blocked

      public double blocked
      Total blocked mass; station queues sum to N minus this.
    • drain

      public Matrix drain
      Rate each waiting queue drains at, one per binding cap.
    • capacityStaged

      public boolean[] capacityStaged
      True where a cap stages the blocked job rather than holding or losing it.
    • capacityRegion

      public int[] capacityRegion
      Which region each cap came from, -1 for a station buffer.
    • capacityStation

      public int[] capacityStation
      Which station each cap came from, -1 for a region cap.
    • capacitySwitches

      public List<double[]> capacitySwitches
      Every time the trajectory made a cap start or stop binding, as {t, row, kind} with kind 0 = release, 1 = activate, 2 = a crossing the cap could not hold. Empty for a steady-state solve.
  • Constructor Details

    • DaeAnalyzer

      public DaeAnalyzer()
  • Method Details