Delay Differential-Algebraic Equations (DDAEs)

• DDAEs are systems that couple differential equations, algebraic constraints, and delayed arguments to model complex dynamic phenomena.
• They are analyzed through spectral, index-based, and geometric methods to ensure solvability, stability, and admissible histories.
• Applications span control design, hybrid numerical-experimental systems, and biological dynamics, driving advanced numerical and structural algorithms.

Delay differential-algebraic equations (DDAEs) are systems in which differential dynamics, algebraic constraints, and delayed arguments are coupled in a single model. Representative forms include the linear time-invariant equation

$E\dot{x}(t)=Ax(t)+Bx(t-\tau)+f(t),$

the nonlinear initial-trajectory problem

$$0=F\bigl(t,z(t),\dot z(t),z(t-\tau)\bigr), \qquad z(t)=\phi(t)\ \text{for}\ t\in[-\tau,0],$$

and algebraic-delay systems

$x'(t)=G(r(t),x_t), \qquad 0=\Delta(r(t),x_t),$

in which the delay is determined only implicitly by an extra algebraic equation. The class also includes systems with state-dependent delay governed by an algebraic equation, systems with distributed delay, and DAEs with infinite delay. Across these settings, the central themes are solvability, admissibility of history functions, spectral and index-based structure, discontinuity propagation, geometric solution manifolds, and numerical regularization (Ha, 2018, Unger, 2020, Walther, 2023, Bisconti et al., 2011).

1. Canonical models and subclasses

A basic linear DDAE is written as

$E\dot{x}(t)=Ax(t)+Bx(t-\tau)+f(t),$

with $E,A,B\in\mathbb{R}^{\ell\times n}$ and $\tau>0$. In the control-oriented literature, interconnected delayed systems are also written in the form

$$\left\{ \begin{array}{l} E\, \dot x(t) = A_0 x(t) + \sum_{i=1}^m A_i x(t-\tau_i) + B w(t),\ z(t) = C x(t), \end{array} \right.$$

where $E$ may be singular, so algebraic constraints are retained rather than eliminated. This representation is used to encode delays in states, inputs, and outputs, as well as direct feedthroughs and controller delays (Ha, 2018, Gumussoy et al., 2020, Gumussoy et al., 2020).

Nonlinear DDAEs arise naturally in hybrid numerical-experimental systems and are written as

$$0 = F\bigl(t, z(t), \dot{z}(t), z(t-\tau)\bigr), \quad t > 0,$$

with a prescribed initial trajectory $z(t)=\phi(t)$ on $$0=F\bigl(t,z(t),\dot z(t),z(t-\tau)\bigr), \qquad z(t)=\phi(t)\ \text{for}\ t\in[-\tau,0],$$0. In real-time dynamic substructuring, a numerical subsystem and an experimental subsystem are each modeled as nonlinear DAEs and then coupled through a transfer system that induces a time-delay $$0=F\bigl(t,z(t),\dot z(t),z(t-\tau)\bigr), \qquad z(t)=\phi(t)\ \text{for}\ t\in[-\tau,0],$$1, producing a nonlinear DDAE after interconnection (Unger, 2020).

A distinct subclass consists of state-dependent delays governed by algebraic equations. One model has the form

$$0=F\bigl(t,z(t),\dot z(t),z(t-\tau)\bigr), \qquad z(t)=\phi(t)\ \text{for}\ t\in[-\tau,0],$$2

where the delay $$0=F\bigl(t,z(t),\dot z(t),z(t-\tau)\bigr), \qquad z(t)=\phi(t)\ \text{for}\ t\in[-\tau,0],$$3 is implicit and nonlocal. Closely related algebraic-delay systems are studied through

$$0=F\bigl(t,z(t),\dot z(t),z(t-\tau)\bigr), \qquad z(t)=\phi(t)\ \text{for}\ t\in[-\tau,0],$$4

with $$0=F\bigl(t,z(t),\dot z(t),z(t-\tau)\bigr), \qquad z(t)=\phi(t)\ \text{for}\ t\in[-\tau,0],$$5 an auxiliary function determined by the algebraic constraint (Hu, 2017, Walther, 2023).

The scope of the field extends further to Volterra delay-integro-differential-algebraic equations (DIDAEs), which add integral operators, and to perturbed DAEs with distributed and possibly infinite delay,

$$0=F\bigl(t,z(t),\dot z(t),z(t-\tau)\bigr), \qquad z(t)=\phi(t)\ \text{for}\ t\in[-\tau,0],$$6

where $$0=F\bigl(t,z(t),\dot z(t),z(t-\tau)\bigr), \qquad z(t)=\phi(t)\ \text{for}\ t\in[-\tau,0],$$7 for $$0=F\bigl(t,z(t),\dot z(t),z(t-\tau)\bigr), \qquad z(t)=\phi(t)\ \text{for}\ t\in[-\tau,0],$$8 and $$0=F\bigl(t,z(t),\dot z(t),z(t-\tau)\bigr), \qquad z(t)=\phi(t)\ \text{for}\ t\in[-\tau,0],$$9 acts on the entire past trajectory in the phase space $x'(t)=G(r(t),x_t), \qquad 0=\Delta(r(t),x_t),$0 (Wang et al., 1 Apr 2025, Bisconti et al., 2011).

2. Solvability, regularity, and admissible histories

For linear time-invariant DDAEs, solvability is characterized spectrally through the regularity of the matrix triple $x'(t)=G(r(t),x_t), \qquad 0=\Delta(r(t),x_t),$1. The characteristic polynomial is

$x'(t)=G(r(t),x_t), \qquad 0=\Delta(r(t),x_t),$2

and the triple is regular if $x'(t)=G(r(t),x_t), \qquad 0=\Delta(r(t),x_t),$3 is not identically zero. With a consistent and sufficiently smooth initial function, the DDAE admits a unique piecewise solution if and only if the triple $x'(t)=G(r(t),x_t), \qquad 0=\Delta(r(t),x_t),$4 is regular. In this setting, a solution means that $x'(t)=G(r(t),x_t), \qquad 0=\Delta(r(t),x_t),$5 is piecewise continuously differentiable, $x'(t)=G(r(t),x_t), \qquad 0=\Delta(r(t),x_t),$6 is continuous, and the equation holds for all $x'(t)=G(r(t),x_t), \qquad 0=\Delta(r(t),x_t),$7 (Ha, 2018).

For nonlinear DDAEs, solvability is tied to index structure and to admissibility of the history. In the strangeness-index framework, the initial trajectory problem

$x'(t)=G(r(t),x_t), \qquad 0=\Delta(r(t),x_t),$8

has local existence and uniqueness if the DDAE is sufficiently smooth, has strangeness index $x'(t)=G(r(t),x_t), \qquad 0=\Delta(r(t),x_t),$9, satisfies Hypothesis 2 with characteristics $E\dot{x}(t)=Ax(t)+Bx(t-\tau)+f(t),$0 and also for $E\dot{x}(t)=Ax(t)+Bx(t-\tau)+f(t),$1, is not advanced, and the initial trajectory $E\dot{x}(t)=Ax(t)+Bx(t-\tau)+f(t),$2 is admissible. Admissibility is expressed by the explicit condition

$E\dot{x}(t)=Ax(t)+Bx(t-\tau)+f(t),$3

This makes the regularity of the history an intrinsic part of the model, not a peripheral compatibility condition (Unger, 2020).

In hybrid numerical-experimental systems, the direct delayed coupling often fails standard regularity conditions. A key structural observation is that shifting certain equations in time can restore regularity and thus solvability. For linear subsystems, the hybrid system has shift index one. When both numerical and experimental subsystems are strangeness-free, the shifted hybrid DDAE also has strangeness index $E\dot{x}(t)=Ax(t)+Bx(t-\tau)+f(t),$4; when the subsystems have higher index, the differentiation index of the hybrid system is at most the sum of the subsystem indices plus one. The same analysis emphasizes that only non-advanced DDAEs, that is retarded or neutral systems, admit well-posed initial trajectory problems in general (Unger, 2020).

3. Spectral structure, stability, and discontinuity propagation

A recurrent misconception is that the eigenvalue-based stability criterion for ordinary DDEs extends unchanged to DDAEs. For linear homogeneous DDAEs, the literature shows that this is only valid for non-advanced systems. In particular, the null solution of

$E\dot{x}(t)=Ax(t)+Bx(t-\tau)+f(t),$5

is exponentially stable if and only if the DDAE is non-advanced and its spectrum lies in the open left half-plane. Counterexamples show that a DDAE can have spectrum in $E\dot{x}(t)=Ax(t)+Bx(t-\tau)+f(t),$6 and still fail to be stable in the $E\dot{x}(t)=Ax(t)+Bx(t-\tau)+f(t),$7-sense because delayed derivatives or high-index effects are present (Ha, 2018).

To account for this, weak exponential stability is introduced. The null solution is $E\dot{x}(t)=Ax(t)+Bx(t-\tau)+f(t),$8-weakly exponentially stable if there exist $E\dot{x}(t)=Ax(t)+Bx(t-\tau)+f(t),$9 and constants $E,A,B\in\mathbb{R}^{\ell\times n}$0 such that

$E,A,B\in\mathbb{R}^{\ell\times n}$1

For DDAEs with pairwise-commuting coefficients, spectral conditions imply $E,A,B\in\mathbb{R}^{\ell\times n}$2-weak exponential stability, where $E,A,B\in\mathbb{R}^{\ell\times n}$3 is the nilpotency index of the relevant nilpotent block. This separates decay properties from the regularity class in which decay can be asserted (Ha, 2018).

A complementary viewpoint classifies DDAEs by the propagation of primary discontinuities under the method of steps. Using the quasi-Weierstraß form for regular matrix pencils, three propagation types are distinguished.

Propagation type	Matrix criterion	Consequence
Smoothing	$E,A,B\in\mathbb{R}^{\ell\times n}$4, $E,A,B\in\mathbb{R}^{\ell\times n}$5 nilpotent	regularity improves after delays
Discontinuity invariant	$E,A,B\in\mathbb{R}^{\ell\times n}$6, $E,A,B\in\mathbb{R}^{\ell\times n}$7 not nilpotent	derivative discontinuity persists
De-smoothing	$E,A,B\in\mathbb{R}^{\ell\times n}$8	regularity can deteriorate

This classification covers all regular DDAEs, takes arbitrary inhomogeneities and history functions as a worst-case baseline, and reveals possible hidden delays in the DDAE. It also refines earlier classifications: the new de-smoothing class is a subset of the previously termed advanced class, so some systems previously grouped as advanced can be structurally benign with respect to discontinuity growth. Splicing conditions, that is matching derivatives of history and solution at transition points, can alter the observed propagation behavior and in some index-3 cases support unique solvability under higher-order matching (Unger, 2018).

Spectral criteria also appear in specialized DDAE kernel problems. For the delay differential equation governing the Lyapunov matrix kernel $E,A,B\in\mathbb{R}^{\ell\times n}$9 of a linear time-delay system with constant and distributed delay, existence and uniqueness are characterized exactly by the condition that the spectrum does not contain opposite pairs: $\tau>0$0 The criterion is necessary and sufficient for the associated auxiliary boundary-value system and hence for the original DDE with integral initial constraints, which the paper interprets as a DDAE setting (Gumussoy et al., 2018).

4. Analytic and geometric reformulations

One analytic strand of DDAE theory replaces delay-algebraic structure by a finite-dimensional boundary-value problem. In the construction of quadratic Lyapunov functionals for a linear time-delay system with a constant and a distributed delay,

$\tau>0$1

the kernel $\tau>0$2 satisfies a DDE with symmetry and algebraic coupling conditions. The problem is reformulated as a delay-free auxiliary ODE system for matrices $\tau>0$3, subject to algebraically coupled split-boundary conditions at $\tau>0$4 and $\tau>0$5. The solution is reconstructed by

$\tau>0$6

The ODE system is linear, finite-dimensional with $\tau>0$7, and time-invariant. In Kronecker form,

$\tau>0$8

and the matrix $\tau>0$9 is invertible exactly when the spectral condition holds. This yields a tractable analytic solution framework for a DDE arising from a DDAE kernel problem (Gumussoy et al., 2018).

A geometric strand studies compatible initial data as a solution manifold. For algebraic-delay systems

$$\left\{ \begin{array}{l} E\, \dot x(t) = A_0 x(t) + \sum_{i=1}^m A_i x(t-\tau_i) + B w(t),\ z(t) = C x(t), \end{array} \right.$$0

the solution manifold is

$$\left\{ \begin{array}{l} E\, \dot x(t) = A_0 x(t) + \sum_{i=1}^m A_i x(t-\tau_i) + B w(t),\ z(t) = C x(t), \end{array} \right.$$1

Under suitable smoothness assumptions, this set is a $$\left\{ \begin{array}{l} E\, \dot x(t) = A_0 x(t) + \sum_{i=1}^m A_i x(t-\tau_i) + B w(t),\ z(t) = C x(t), \end{array} \right.$$2 submanifold of codimension $$\left\{ \begin{array}{l} E\, \dot x(t) = A_0 x(t) + \sum_{i=1}^m A_i x(t-\tau_i) + B w(t),\ z(t) = C x(t), \end{array} \right.$$3 in $$\left\{ \begin{array}{l} E\, \dot x(t) = A_0 x(t) + \sum_{i=1}^m A_i x(t-\tau_i) + B w(t),\ z(t) = C x(t), \end{array} \right.$$4. Hans-Otto Walther shows that there exists a diffeomorphism

$$\left\{ \begin{array}{l} E\, \dot x(t) = A_0 x(t) + \sum_{i=1}^m A_i x(t-\tau_i) + B w(t),\ z(t) = C x(t), \end{array} \right.$$5

that straightens the manifold to a simpler form on a strip $$\left\{ \begin{array}{l} E\, \dot x(t) = A_0 x(t) + \sum_{i=1}^m A_i x(t-\tau_i) + B w(t),\ z(t) = C x(t), \end{array} \right.$$6, where

$$\left\{ \begin{array}{l} E\, \dot x(t) = A_0 x(t) + \sum_{i=1}^m A_i x(t-\tau_i) + B w(t),\ z(t) = C x(t), \end{array} \right.$$7

For explicit delay systems, the transformed manifold is a graph. The result extends previous work by removing the restriction that the implicitly defined delay must be strictly negative, so delays may vanish (Walther, 2023).

For infinite delay, Bisconti and Spadini analyze

$$\left\{ \begin{array}{l} E\, \dot x(t) = A_0 x(t) + \sum_{i=1}^m A_i x(t-\tau_i) + B w(t),\ z(t) = C x(t), \end{array} \right.$$8

with $$\left\{ \begin{array}{l} E\, \dot x(t) = A_0 x(t) + \sum_{i=1}^m A_i x(t-\tau_i) + B w(t),\ z(t) = C x(t), \end{array} \right.$$9 acting on the entire past history. Under the kernel conditions

$E$0

and an invertibility condition on the algebraic block of $E$1, the system is equivalent to a retarded functional differential equation on a manifold. If the Brouwer degree of the vector field is well defined and nonzero on an appropriate open set, there exists a connected set $E$2 of non-trivial $E$3-periodic solution pairs $E$4, and $E$5 is not compact in the set of trivial solutions. This gives a global continuation result for periodic solutions in a DAE with infinite delay (Bisconti et al., 2011).

5. Structural algorithms and numerical methods

Structural preprocessing for DDAEs generalizes the Pantelides algorithm from DAEs to systems with delay. The central insight is that differentiation and shifting are very similar from a structural point of view. In a graph-theoretical formulation, equations and variables are represented in a bipartite graph, while variables with different derivative order or time shift are linked through equivalence classes. This allows the algorithm to decide which equations need to be differentiated or shifted in order to construct a solution. A necessary and sufficient criterion for termination is obtained: the Pantelides DDAE algorithm terminates if and only if the system is structurally nonsingular with respect to the equivalence classes $E$6. The framework also includes differentiation during shifting, trimmed linearization, and update rules for the graph as operations are performed (Ahrens et al., 2019).

For nonlinear Volterra DIDAEs with constant delay, numerical stability has been studied through hybrid schemes that combine Runge-Kutta methods with compound quadrature rules. The resulting CQRK methods support algebraic constraints and are analyzed in stiff regimes, including cases with large classical Lipschitz constants. The exact solution is treated through a generalized Halanay inequality, yielding exponential decay under

$E$7

For the numerical method, global stability is defined by bounded dependence on the initial perturbation, while asymptotic stability requires

$E$8

Theorems 1 and 2 then give constructive conditions involving algebraic stability of the underlying Runge-Kutta method, a root condition $E$9, the step size $$0 = F\bigl(t, z(t), \dot{z}(t), z(t-\tau)\bigr), \quad t > 0,$$0, the delay $$0 = F\bigl(t, z(t), \dot{z}(t), z(t-\tau)\bigr), \quad t > 0,$$1, quadrature constants, and the Lipschitz or one-sided Lipschitz constants. Numerical examples use the 2-stage Lobatto III C Runge-Kutta method with Simpson's quadrature and show perturbation decay consistent with the theory (Wang et al., 1 Apr 2025).

These developments indicate that DDAE numerics depends on two coupled tasks: structural regularization of the delayed-algebraic model and stability-preserving discretization once a suitable form has been obtained. The literature does not treat these tasks as interchangeable; rather, the structural graph analysis and the numerical stability analysis are complementary stages.

6. Applications in hybrid systems, control, and bifurcation

In real-time dynamic substructuring, a complex structure is subdivided into a numerical part and an experimental part, coupled in real time through a transfer system. The transfer system introduces a time-delay $$0 = F\bigl(t, z(t), \dot{z}(t), z(t-\tau)\bigr), \quad t > 0,$$2, and the resulting hybrid model is typically described by nonlinear DDAEs. The analysis shows how shifting restores regularity, how admissible initial trajectories must satisfy delayed consistency conditions, and why the distinction between retarded, neutral, and advanced structure matters for practical well-posedness. The motivating example is a coupled pendulum–mass–spring–damper system, and the paper notes that port-Hamiltonian subsystems are particularly advantageous because their structure ensures good index and robust solvability properties (Unger, 2020).

In robust control, DDAEs provide a systematic model class for plants and controllers with delays in state, input, and output, without using elimination techniques. For the transfer function

$$0 = F\bigl(t, z(t), \dot{z}(t), z(t-\tau)\bigr), \quad t > 0,$$3

the standard $$0 = F\bigl(t, z(t), \dot{z}(t), z(t-\tau)\bigr), \quad t > 0,$$4 norm may be sensitive to arbitrarily small delay perturbations. To address this, the strong $$0 = F\bigl(t, z(t), \dot{z}(t), z(t-\tau)\bigr), \quad t > 0,$$5 norm is introduced as

$$0 = F\bigl(t, z(t), \dot{z}(t), z(t-\tau)\bigr), \quad t > 0,$$6

with the key identity

$$0 = F\bigl(t, z(t), \dot{z}(t), z(t-\tau)\bigr), \quad t > 0,$$7

A predictor-correct algorithm computes the norm, and gradients of the strong norm with respect to controller parameters support non-smooth, non-convex optimization for fixed-order or structured controller design. The associated analysis concludes that the strong $$0 = F\bigl(t, z(t), \dot{z}(t), z(t-\tau)\bigr), \quad t > 0,$$8 norm is the more appropriate performance metric in practical control applications whenever high-frequency paths are present in control loops (Gumussoy et al., 2020, Gumussoy et al., 2020).

Bifurcation theory for DDAEs with implicit state-dependent delay uses the $$0 = F\bigl(t, z(t), \dot{z}(t), z(t-\tau)\bigr), \quad t > 0,$$9-equivariant degree. For

$z(t)=\phi(t)$0

a local Hopf bifurcation occurs when an isolated center has nonzero crossing number, and a global Hopf theorem states that the connected component of nontrivial periodic solutions is either unbounded in the product of state, parameter, and period space, or meets only finitely many other bifurcation points, with the sum of signed crossing numbers vanishing. The theory is applied to an extended Goodwin model in which the threshold-type algebraic delay

$z(t)=\phi(t)$1

may vanish at the stationary state, a case for which standard reductions do not apply. This suggests that DDAEs are not only a unifying modeling framework but also a natural setting for global continuation results in implicitly delayed biological dynamics (Hu, 2017).