Abstract
The focus of this paper is on the use of polynomial chaos (PC) for developing surrogate models for differential algebraic equations (DAEs) with time-invariant uncertainties. Intrusive and nonintrusive approaches to synthesize PC surrogate models are presented including the use of Lagrange interpolation polynomials as basis functions. Unlike ordinary differential equations (ODEs), if the algebraic constraints are a function of the stochastic variable, some initial conditions of the DAEs are also random. A benchmark RLC circuit which is used as a benchmark for linear models is used to illustrate the development of a PC-based surrogate model. A nonlinear example of a simple pendulum also serves as a benchmark to illustrate the potential of the proposed approach. Statistics of the results of the PC models are validated using Monte Carlo (MC) simulations in addition to estimating the evolving probably density functions (PDFs) of the states of the pendulum.
1 Introduction
A variety of dynamical systems are governed by a mixture of differential and algebraic equations referred to as differential algebraic equations (DAEs). Musau et al. [1] list and describe benchmark linear DAE problems including mechanical systems, electrical RLC circuits, and electrical generator system, among others. Campbell et al. [2] focus on optimal control design for DAE problems including tracking problems, path constraint problems, servo-constraints approach in mechanical systems, and other engineering problems. They discuss formal methods of index reduction of the DAEs and numerical algorithms available for solving these practical DAE problems. Similar to ordinary differential equations (ODEs), identifying analytical solutions for all DAEs are seldom possible. Hence, robust numerical techniques are required to compute solutions for DAEs. A general approach for finding a solution for a DAE is to reduce its index at the onset (a detailed discussion on indices of DAEs is provided later). Next a classification, whether a DAE is a boundary value problem (BVP) or an initial value problem (IVP) is performed. Depending on the type of DAE (BVP or IVP), different numerical algorithms are used to compute solution for the DAE problem. Whenever a numerical method is proposed for solving a mathematical problem, it is necessary to discuss the stability, convergence, and other standard metrics for the respective numerical method. Lamour et al. [3] have compiled a survey on BVPs for DAEs. They first provide the analytical theory for BVPs in the DAEs' framework. Later, they present collocation methods for well-posed BVPs, shooting method, and other miscellaneous numerical methods for solving BVPs in DAEs. Brenan et al. [4] have discussed basic theory for DAEs, solvability, index, and index reduction method for both linear and nonlinear DAEs. Furthermore, they have discussed multistep and one-step methods for solving DAEs. A software for solving DAEs, dassl is discussed. Ilchmann and Reis [5] have also extensively studied DAEs. Their paper starts with a detailed discussion of the general theory of linear DAEs and their observability. Model reduction of chemical processes and several optimal control problems posed in a DAE framework are discussed. Further, a functional-analytic overview of DAE systems is presented. Ascher and Petzold [6] have also given an overview on general mathematical structure for DAEs. Moreover, numerical methods (backward Euler methods, backward differentiation formula and multistep method, Radau collocation and implicit Runge–Kutta methods, and specialized Runge–Kutta method for Hessenberg index two DAEs) for solving DAEs are proposed, and stability, convergence, error estimation, and other standard metric of numerical methods are discussed. Brown et al. [7] proposed a new algorithm for the solution of large-scale systems of DAEs. Unlike IVPs for ODEs, initial conditions for the system states of DAEs require special attention and treatment. Brown et al. [8] have described a new algorithm for the evaluation of consistent initial conditions for a class of DAEs. Based on minimum information provided by the user, the algorithm calculates consistent initial conditions, which are required at the beginning of any numerical integration technique for solving DAEs. Hence, there is considerable literature available on the development of both mathematical theorem and numerical solution algorithms for DAE systems [9,10].
Polynomial chaos (PC) expansion is a well-established methods for uncertainty quantification (UQ) [11–16]. There are two major approaches to synthesizing PC surrogate models for systems with uncertainties, intrusive PC (example: Galerkin projection) and nonintrusive PC (example: collocation method) method. While intrusive method requires analytical models of the uncertain system, nonintrusive method can circumvent the necessity of such information. Hence, nonintrusive PC methods draw attention in situations where the model can be considered a black-box. Kim et al. [12] have presented a review of both intrusive and nonintrusive formulation. Eldred and Burkardt [13] provide a comprehensive review of the mathematical formulation and application of nonintrusive PC.
Pulch [17–21] has shown in a very well structured way how PC technique (intrusive Galerkin projection method and nonintrusive collocation method) can be used for UQ of a variety of stochastic DAE systems. He has presented both the collocation and intrusive Galerkin projection approaches to quantify uncertainty in multirate partial DAEs [17]. PC yields a larger number of deterministic equations from the original stochastic equation. Pulch has presented a systematic way of calculating the index of the resulting deterministic DAEs [19,20]. Further, he presents a generalized approach for using Galerkin projection-based PC for analyzing periodic processes of stochastic DAEs [18]. In this paper, in contrast to the Wiener–Askey scheme, Lagrange interpolation polynomials are used as basis functions. This novel approach to parameterizing the polynomial chaos expansion results in a sampling-based approach for the development of the PC surrogate model, which are identical to those developed via the intrusive approach for linear systems.
Consistent evaluation of initial conditions and boundary conditions is imperative in obtaining accurate solutions for a DAE system. With the increase in the number of equations yielded from PC expansion, evaluation of consistent initial conditions becomes complex. In this work, the authors have attempted to address this problem. In this work, with the help of benchmark examples of stochastic DAEs, the authors have shown concepts like index reduction of higher-order DAEs, identification of hidden constraints, and evaluation of consistent initial conditions.
The paper is organized as follows: Sec. 2 of the paper provides a general overview of DAEs to assist the reader in comprehending the basic concepts and forms of DAEs. In Sec. 3, a general overview of UQ for a system of DAEs using PC expansion is discussed. Finally, in Sec. 4, benchmark DAE examples are solved using both intrusive and nonintrusive PC method, and results illustrating the ability of PC models to capture the statistics of the evolving uncertain states are presented.
2 General Overview of Differential Algebraic Equations
where the Jacobian of G with respect to (denoted by ) is nonsingular, or this may also occur because in Eq. (1) is singular. These systems are referred to as differential algebraic equations.
where we can distinguish between the differential variables x(t) and the algebraic variables z(t). Note that Eqs. (4a), (4b), (5a), and (5b) are both semi-explicit DAEs.
By repeated differentiation of Eq. (5b) and by substitution, one can yield an explicit set of ODEs. The number of derivatives required before the explicit set of ODEs are realized is called the index of the DAEs. (Note, ODEs have an index of 0.) The concept of index is very important characteristics for understanding the complexity of the system of DAEs. In general, higher the index of the DAEs, harder it is to find the solution of the DAEs.
To reiterate, in case of the semi-explicit form of DAEs, the index of the DAEs is the minimum number of times the constraint Eq. (4b) needs to be differentiated (let's say times) so that the and together form a system of ODEs which is required to find an unique solution for . From the study of ODEs, it can be said that the initial condition for a vector valued is independent of each other. In contrast, the initial conditions of the state variables for DAEs have to be consistent. (Note that in some problems, the states are not all prescribed at the same value of the independent variables and can result in boundary conditions rather than initial conditions.) In other words, the initial conditions and boundary values must satisfy the algebraic constraints of the DAEs for all values of the independent variable. In index one semi-explicit DAEs, the initial conditions only satisfy the explicitly stated algebraic constraint. In semi-explicit DAEs with higher-order indices, in addition to the explicit algebraic constraint (4b) and (5b), hidden constraints need to be enforced. Every equation of the set of equations derived after differentiation () of is called a hidden constraint.
Equation (6b) yields two solution for : and .
Setting , Eq. (6c) yields . After substituting the solution for in Eq. (6a), the solution for yields (t0 is initial time). The system of DAEs has index one, and the complete solution is Hence, the system of DAEs is an index one system in semi-explicit form.
Setting , Eq. (6c) yields . After differentiating Eq. (6c) once, the first equation results in . The system of DAEs has index two, and the complete solution is . Hence, the system of DAEs is an index two system in semi-explicit form.
The solution to Eq. (7) is and . It needs to be emphasized that both the solution should still satisfy Eq. (6b). However, does not satisfy Eq. (6b), which is an algebraic constraint in the system of DAE, Eqs. (6a)–(6c). Hence, the only feasible solution is .
where is the Jacobian of with respect to the vector x. If is nonsingular, the system (Eq. (11)) is an implicit ODE, and it can be said that Eq. (10) has index one. If is singular, that implies, Eq. (11) continues to include algebraic constraints and additional derivatives are required before a set of ODEs' result. For instance, if after evaluating the derivative of Eq. (11), the coefficient matrix of the highest derivative of y(t) is nonsingular, the original problem has index two. One can define higher indices similarly.
Most of the higher index DAE problems observed in nature can be expressed as a combination of a system of coupled and nonlinear ODEs and a system of coupled and nonlinear algebraic constraints. In such DAEs, the classification of algebraic state variables and differential variables might not be explicitly distinguishable. Moreover, for some higher index DAEs of special forms, the algebraic variables may all be eliminated (in principle) using the same number of differentiations as its index. These DAEs of the special forms are called forms of the DAEs as described in Ref. [6].
In this paper, the RLC benchmark problem described in Ref. [22] has been chosen as an example of linear DAEs with constant coefficients, and a simple pendulum is chosen as an example of nonlinear DAEs to illustrate the proposed uncertainty quantification approach.
3 Uncertainty Quantification Using Polynomial Chaos
The aforementioned development catered to deterministic DAEs. Often the parameters which define the models are uncertain and are characterized by probability distributions. There is consequently a need to develop tools for the characterization of the statistics of the evolving states of stochastic DAEs. Polynomial chaos is a well-established approach for the development of surrogate models for ODEs which provides a computationally efficient approach to estimate the mean, variance, and higher moments of the evolving states of the stochastic DAE models.
3.1 General Overview of Using Polynomial Chaos in Differential Algebraic Equations Framework.
will be developed, where is a random vector with probably density function (PDF) . Xiu and Karniadakis [23] have presented a list of orthogonal polynomial to use in the PC expansion associated with the most common distribution functions. It is commonly known as scheme. Depending on the uncertainties present in the system parameters, a vector of random variables with different probability distribution functions can exist. For the sake of simplicity in explaining UQ (using PC), only one random variable, , is chosen and the method is explained.
For computational tractability, the PC expansion is truncated to a finite order which results in an approximation of the exact solution.
Hence, for N number of orthogonal bases, the number of initial equations () become number of equations ( stands for multiplication operation here), which are needed to be solved simultaneously in order to find number of unknown functions of time.
Higher moments can be similarly calculated and are studied in detail by Lefebvre [24].
3.2 Nonintrusive Polynomial Chaos With Lagrange Polynomial Bases.
Following the conception of homogeneous chaos by Wiener for Gaussian processes, the idea was generalized by Xiu and Karniadakis using the Wiener–Askey scheme to cater to a large class of distributions. Nandi and Singh [25] illustrated that using Lagrange polynomial functions as basis function to parameterize the PC expansion for linear systems results in a sampling-based approach to develop a surrogate model. They show that the resulting surrogate model is identical to that resulting from an intrusive Galerkin projection-based approach. This precludes the need for multivariate integrals making it an attractive approach for developing PC surrogate models where the system of interest can be considered a black-box.
where λi is the Christoffel number, and the quadrature points ζi are the roots of the Nth order polynomial orthogonal with respect to the distribution . For example, for a uniformly distributed variable, ζi are the roots of the Legendre polynomial.
is exact. The quadrature points ζj are the roots of the Legendre polynomial of degree n. Polynomial chaos expansion using Lagrange polynomials as basis functions will be used to develop surrogate models for the two examples presented in Sec. 4.
4 Numerical Examples
Two benchmark problem will be considered in this section to illustrate the use of PC to characterize the evolution of uncertainties of the states of the system. The first model is an RLC circuit which is a linear model, and the second problem is a simple pendulum modeled in Cartesian coordinates which is a nonlinear problem.
4.1 Benchmark Problem of a Single Loop RLC Circuit.
In these DAEs, there are four system states (). In general, ( denotes as equivalent to) is a function of time and is a known input.
variable is also an index one algebraic variable which can also be eliminated simply by substitution. The system described by Eq. (27) will be referred to as in this work.
After solving the system of the ODEs, and can be computed from Eqs. (26c) and (27c), respectively. The system described by Eq. (28) will be referred to as in this work.
The system described by Eq. (29) will be referred to as in this work. does not require any enforcement of algebraic constraint given by Eq. (27c) at the initial time. Initial conditions for both states ( and ) can be mentioned as independent initial conditions. On the contrary, does require enforcement of aforementioned algebraic constraint at the initial time. Any two of the three states can be assigned independent initial conditions, while the third state will have an initial condition which is consistent with (Eq. (27c)).
4.1.1 Enforcing Initial Condition Constraint.
and are Legendre polynomials, and denotes the dot product of two vectors.
which are 2N deterministic equations which are a function of the deterministic initial conditions and .
which results in 3N deterministic equations which are function of the initial conditions given in Eqs. (32a)–(32c).
4.1.2 Numerical Results
System model.
The three models , and will be studied with , and . The initial conditions in conjunction with the uncertain resistance make an uncertain variable. The coefficients of the PC expansion of are determined via Eq. (32c). The DAE is solved using the “ode15i” function in matlab. To validate the results of the PC model, 10,000 Monte Carlo (MC) simulations are used to estimate the mean and variance of the evolving uncertain states and will be considered the “truth.” For the PC expansion, a sixth-order expansion is used to generate the following results.
Figure 2(a) illustrates the time evolution of the mean of the variables I(t), , and which are given by the first coefficient of the polynomial chaos expansion of the same variables. It is clear that the mean estimated via the PC and the MC is coincident. We arrive at the same conclusion by examining the time evolution of the variances (Fig. 2(b)). It should be noted that the variance of is nonzero at the initial time reflecting the uncertainty in the initial value of the current generated by the uncertainty in the resistance.

Comparison of PC and Monte Carlo results from model. (a) Estimate of state means and (b) estimate of state variances.
System model.
In the attempt to reduce the index of the DAE, a simple method of substitution was used (as explained in Eqs. (28a) and (28b)) but it should be kept in mind that, although the substitution of algebraic variables in the DAE may appear intuitive in case of linear DAE, the procedure can be complex in many cases when handling nonlinear DAEs. The system model for is solved by using “” function of matlab. Sixth-order PC expansion was used to represent the system states and .
This model implicitly satisfies the algebraic constraint on the system states at every instances of time. The results in the figures show that the means and variances of the system states calculated using PC formulation exactly match the ones derived using MC formulation. Since there is a constraint on the initial conditions of the system states, nonzero variance (Fig. 3(b)) for at the initial time was observed.

Comparison of PC and Monte Carlo results from model. (a) Estimate of state means and (b) estimate of state variances.
System model.
In the attempt to reduce the index of the DAE (the method explained in Eqs. (29a)–(29c) in detail), the differentiation of the algebraic equation followed by appropriate substitution and rearrangement was implemented. A sixth-order PC expansion was used to develop a surrogate model for the stochastic ODEs, and Eqs. (37a)–(37c) were solved in order to compute the time-dependent PC coefficients of the respective system states. The calculation for the dependent initial condition, i.e., , is provided in Eq. (35). The ODEs were solved using ode45 function of matlab. The number of basis selection for PC expansion and sampling size for MC calculation remains same as the DAE's formulation.
It can be observed from Fig. 4 that the means and variances estimated from the PC model and MC simulation are coincident and match the results from the previous simulations of systems and .

Comparison of PC and Monte Carlo results from model. (a) Estimate of state means and (b) estimate of state variances.
4.2 Simple Pendulum.
It can be observed that the DAEs are not of semi-explicit form. The algebraic constraint in Eq. (48b) needs to be differentiated two times in order to express it into a fully implicit ODE.
for all the instances of time t.
since A is a nonsingular matrix for every instance of time t.
4.2.1 Treatment of the Initial Conditions.
Unlike Eqs. (37a)–(37c) where Galerkin projection resulted in closed-form expressions, the same is not possible for Eqs. (69a) and (69b) which will be referred to as . Consequently, a nonintrusive PC technique is implemented to characterize the statistics of the stochastic states. The specific nonintrusive approach presented here is referred to as the stochastic response surface method which uses a linear least squares approach to solve for the coefficients of the PC expansion. A set of response values is obtained for specific realizations of the stochastic parameters, where more samples (M) than the number of PC coefficients (N) to be determined are used to pose an overdetermined least squares problem.
Note that the pseudo-inverse term needs to only be evaluated once. A similar equation is derived to determine the PC coefficients of the .
The Gauss quadrature formula using N quadrature points ensures an exact evaluation of integrals of polynomials of order . Since the terms in Eqs. (74c) and (74d) have a highest order of , inner product on the N − 1 order basis functions would result in polynomials of order which would only result in the Gauss quadrature formula approximating the inner product integrals.
4.2.2 Numerical Results.
Figure 6 illustrates the time variation of the mean of the four variables listed in Eqs. (76)–(78). It can be noted that the PC results for the DAE and ODE models overlap with the results of the Monte Carlo simulations.
Figure 7 illustrates the time variation of the variance of the four variables listed in Eqs. (76)–(78). It can be noted that the PC results for the DAE and ODE models overlap with the results of the Monte Carlo simulations. It can be noted from the second panel in Fig. 7 that the initial variance of the displacement in y coordinate is nonzero. This reflects the uncertainty in the length of the pendulum which is transferred to the initial condition of the y coordinate when the x coordinate is specified.
To illustrate that the PC surrogate model captures the evolution of the PDF of the evolving states, 1 × 105 number of Monte Carlo simulation of the DAE are used to estimate the PDF of the evolving states using kernel density functions which are illustrated in Fig. 8. Similarly, 1 × 105 number of Monte Carlo samples of the uncertain pendulum length are used to generate samples of the PC states using the surrogate model developed using six basis functions. Results of kernel density estimated PDF of the DAE states are illustrated in Fig. 9. It can be seen that the complex shaped PDFs estimated by the Monte Carlo simulations are captured by the PC surrogate model lending confidence to the use of the PC surrogate model for the DAE. Table 1 presents the Wasserstein distance between the PDFs estimated by the MC and PC approaches which illustrates that the PDFs of all the states at all time are practically identical. Similar evolution of the PDF of the pendulum represented in polar coordinates was completed and shown to match the Monte Carlo results and has not been presented in the interest of brevity.
Wasserstein distance between PDFs estimated from Monte Carlo and polynomial chaos
Wasserstein distance between MC method and surrogate PC model calculated at different instance of time | ||||
---|---|---|---|---|
State | ||||
Time | Z1 | Z2 | Z3 | Z4 |
0.00 s | 9.9565 × 10−7 | 7.1807 × 10−5 | 0 | 0 |
0.45 s | 3.2802 × 10−7 | 5.0317 × 10−5 | 5.7881 × 10−7 | 9.0814 × 10−5 |
0.90 s | 1.0056 × 10−6 | 4.0612 × 10−6 | 3.3365 × 10−6 | 1.1582 × 10−4 |
1.35 s | 1.0185 × 10−6 | 3.4446 × 10−5 | 3.6109 × 10−6 | 9.5807 × 10−5 |
1.80 s | 3.2616 × 10−6 | 7.4939 × 10−5 | 1.583 × 10−5 | 1.4622 × 10−5 |
2.25 s | 7.062 × 10−6 | 5.1427 × 10−5 | 4.7912 × 10−6 | 1.0248 × 10−4 |
2.70 s | 6.9795 × 10−6 | 4.1967 × 10−6 | 5.9642 × 10−5 | 8.0802 × 10−5 |
3.15 s | 2.4678 × 10−5 | 2.2153 × 10−5 | 4.5262 × 10−5 | 3.2469 × 10−5 |
3.60 s | 1.8977 × 10−5 | 4.7335 × 10−5 | 9.7543 × 10−5 | 1.2413 × 10−4 |
4.05 s | 4.5417 × 10−5 | 3.0854 × 10−5 | 2.3585 × 10−5 | 7.7478 × 10−5 |
4.50 s | 3.1358 × 10−5 | 2.4963 × 10−5 | 1.0882 × 10−4 | 3.3921 × 10−4 |
Wasserstein distance between MC method and surrogate PC model calculated at different instance of time | ||||
---|---|---|---|---|
State | ||||
Time | Z1 | Z2 | Z3 | Z4 |
0.00 s | 9.9565 × 10−7 | 7.1807 × 10−5 | 0 | 0 |
0.45 s | 3.2802 × 10−7 | 5.0317 × 10−5 | 5.7881 × 10−7 | 9.0814 × 10−5 |
0.90 s | 1.0056 × 10−6 | 4.0612 × 10−6 | 3.3365 × 10−6 | 1.1582 × 10−4 |
1.35 s | 1.0185 × 10−6 | 3.4446 × 10−5 | 3.6109 × 10−6 | 9.5807 × 10−5 |
1.80 s | 3.2616 × 10−6 | 7.4939 × 10−5 | 1.583 × 10−5 | 1.4622 × 10−5 |
2.25 s | 7.062 × 10−6 | 5.1427 × 10−5 | 4.7912 × 10−6 | 1.0248 × 10−4 |
2.70 s | 6.9795 × 10−6 | 4.1967 × 10−6 | 5.9642 × 10−5 | 8.0802 × 10−5 |
3.15 s | 2.4678 × 10−5 | 2.2153 × 10−5 | 4.5262 × 10−5 | 3.2469 × 10−5 |
3.60 s | 1.8977 × 10−5 | 4.7335 × 10−5 | 9.7543 × 10−5 | 1.2413 × 10−4 |
4.05 s | 4.5417 × 10−5 | 3.0854 × 10−5 | 2.3585 × 10−5 | 7.7478 × 10−5 |
4.50 s | 3.1358 × 10−5 | 2.4963 × 10−5 | 1.0882 × 10−4 | 3.3921 × 10−4 |
5 Conclusions
The focus of this paper is on uncertainty quantification of dynamical system represented as differential algebraic equations. The generalized polynomial chaos approach, which uses orthogonal basis function based on the Wiener–Askey scheme, is used to approximate the stochastic states. Intrusive and nonintrusive methods to develop surrogate models are presented including a nonintrusive approach using Lagrange interpolation polynomials serving as basis functions. A RLC circuit representing a linear system and a simple pendulum represented in Cartesian coordinate which results in a nonlinear equation are used as benchmark problems to study the ability of PC surrogate model to capture the evolving statistics of the uncertain states of the systems. Besides capturing the evolution of the mean and variance, kernel density estimates of the evolving PDFs using Monte Carlo simulations of the DAE and those resulting from the PC surrogate model are used to illustrate the remarkable ability of the PC model to capture the evolving PDFs of the system states. Future work will exploit results of polynomial chaos-based DAEs to solve optimal control problems.
Funding Data
National Science Foundation (NSF) (Award No. CMMI-2021710; Funder ID: 10.13039/100000001).