## 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

*F*and

*y*are vectors. The vector of differential equations is called implicit if the state derivative $y\u2032(t)$ cannot be isolated and represented in terms of the states

*y*. If Eq. (1) can be rewritten in the explicit or canonical form

where the Jacobian of *G* with respect to $y\u2032$ (denoted by $dG/dy\u2032=Gy\u2032$) is nonsingular, or this may also occur because $Fy\u2032$ in Eq. (1) is singular. These systems are referred to as differential algebraic equations.

*a*), (3

*b*), (4

*a*), and (4

*b*). Equations (3

*a*) and (3

*b*) can be rewritten as

where we can distinguish between the *differential variables x*(*t*) and the algebraic variables *z*(*t*). Note that Eqs. (4*a*), (4*b*), (5*a*), and (5*b*) are both semi-explicit DAEs.

By repeated differentiation of Eq. (5*b*) 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. (4*b*) needs to be differentiated (let's say $k$ times) so that the $G(t,y(t))$ and $(dkF(t,y(t)))/dtk$ together form a system of ODEs which is required to find an unique solution for $y(t)$. From the study of ODEs, it can be said that the initial condition for a vector valued $y(t)$ is independent of each other. In contrast, the initial conditions of the state variables $y(t)$ 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 (4*b*) and (5*b*), hidden constraints need to be enforced. Every equation of the set of $k\u22121$ equations derived after $jth$ differentiation ($j=1,2,\u2026k$) of $F(t,y(t))$ is called a hidden constraint.

*a*) and (5

*b*)) but also on the solution of it. It so happens because of the local linearization, i.e., the matrix $Gx(t,x(t),z(t))$ depends on the solution of

*z*(

*t*). In the following example, taken from Chap. 9 of Ref. [6], this observation is explained:

Equation (6*b*) yields two solution for $y2(t)$: $y2(t)=0$ and $y2(t)=1$.

Setting $y2(t)=1$, Eq. (6

*c*) yields $y3(t)=t$. After substituting the solution for $y3(t)$ in Eq. (6*a*), the solution for $y1(t)$ yields $y1(t)=(t2/2)+y1(t0)$ (*t*_{0}is initial time). The system of DAEs has index one, and the complete solution is $y(t)=[y1(t)y2(t)y3(t)]T=[(t2/2)+y1(t0)0t]T$ Hence, the system of DAEs is an index one system in semi-explicit form.Setting $y2(t)=1$, Eq. (6

*c*) yields $y1(t)=t$. After differentiating Eq. (6*c*) once, the first equation results in $y3(t)=1$. The system of DAEs has index two, and the complete solution is $y(t)=[y1(t)y2(t)y3(t)]T=[t11]T$. Hence, the system of DAEs is an index two system in semi-explicit form.

The solution to Eq. (7) is $y2(t)=1/2$ and $y2\u2032(t)=0$. It needs to be emphasized that both the solution should still satisfy Eq. (6*b*). However, $y2(t)=1/2$ does not satisfy Eq. (6*b*), which is an algebraic constraint in the system of DAE, Eqs. (6*a*)–(6*c*). Hence, the only feasible solution is $y2\u2032(t)=0$.

*p*times) the existing DAEs need to be differentiated to solve for $y(t)$ uniquely

*t*, the result is

where $Bx$ is the Jacobian of $B$ with respect to the vector *x*. If $By$ is nonsingular, the system (Eq. (11)) is an implicit ODE, and it can be said that Eq. (10) has index one. If $By$ 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 $Hessenberg$ 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.

*a*) and (3

*b*) and Eqs. (4

*a*) and (4

*b*) can be generalized to represent stochastic DAEs with probabilistic time-invariant model parameters characterizing the system. Equations (3

*a*) and (3

*b*) being the most general form, the polynomial chaos method to quantify uncertainty for that class of DAEs which are represented as

will be developed, where $\zeta $ is a random vector with probably density function (PDF) $P(\zeta )$. 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 $Wiener\u2212Askey$ 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, $\zeta $, 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.

*a*) and (12

*b*) on each of the orthogonal polynomial basis, Eq. (14) result

Hence, for *N* number of orthogonal bases, the number of initial equations ($dim(G)+dim(g)$) become $(dim(G)$$+dim(g))\xd7N$ number of equations ($\xd7$ stands for multiplication operation here), which are needed to be solved simultaneously in order to find $dim(y(t,\zeta ))\xd7N$ 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.

*ζ*is a quadrature point, are orthogonal since they have an inner product

^{i}where *λ _{i}* is the Christoffel number, and the quadrature points

*ζ*are the roots of the

^{i}*N*th order polynomial orthogonal with respect to the distribution $P(\zeta )$. For example, for a uniformly distributed variable,

*ζ*are the roots of the Legendre polynomial.

^{i}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.

for $m=1,2,\u2026,N$. If the terms on the left-hand and right-hand side of Eq. (18) are polynomials of degree less than or equal to $2n\u22121$, then Eq. (24*b*) represents Eq. (24*a*) exactly, else the approximation tends to the exact solution as $n\u2192\u221e$ [26] if the functions in Eq. (18) are continuous.

## 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.

*R*), the inductor (

*L*), and the capacitor (

*C*), respectively. Then, from Kirchoff's laws, the following circuit equations (DAEs) are derived (Fig. 1):

In these DAEs, there are four system states ($I(t),VL(t),VR(t),\u2009and\u2009VC(t)$). In general, $Vs\u225cVs(t)$ ($\u225c$ denotes as *equivalent to*) is a function of time and is a known input.

*c*) yields $VR\u225cVR(I(t),R)$. $VR(t)$ is an index one variable as it can be eliminated simply by substitution (there is no need for differentiation of the variable). After substitution the reduced system of DAEs of index one (detailed discussions of indices of DAEs can be found in Refs. [4] and [6]) are as follows:

variable $VL(t)$ 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 $\mathbb{R}L\u21021$ in this work.

After solving the system of the ODEs, $VR$ and $VL$ can be computed from Eqs. (26*c*) and (27*c*), respectively. The system described by Eq. (28) will be referred to as $\mathbb{R}L\u21022$ in this work.

*a*) and (28

*b*), the system model can be represented as an ODE by differentiating Eq. (27

*c*) and substituting Eqs. (27

*a*) and (27

*b*) in the resulting differential equation, which is represented in terms of three states $I(t),\u2009VC(t)$, and $VL(t)$) as

The system described by Eq. (29) will be referred to as $\mathbb{R}L\u21023$ in this work. $\mathbb{R}L\u21022$ does not require any enforcement of algebraic constraint given by Eq. (27*c*) at the initial time. Initial conditions for both states ($I(t)$ and $VC(t)$) can be mentioned as independent initial conditions. On the contrary, $\mathbb{R}L\u21023$ 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. (27*c*)).

#### 4.1.1 Enforcing Initial Condition Constraint.

and $\varphi i(\zeta )$ are Legendre polynomials, and $\xb0$ denotes the dot product of two vectors.

*c*) onto the random space spanned by the polynomial basis $\varphi k(\zeta )$. This requires triple product integrals with the weighting function $P(\zeta )=0.5$. The identity

*j*symbol is

*c*) to

*a*) and (28

*b*) ($\mathbb{R}L\u21022$) and Eqs. (29

*a*)–(29

*c*) ($\mathbb{R}L\u21023$) will be studied. Galerkin projection of the ODEs of model $\mathbb{R}L\u21022$ after the polynomial expansion of $VC(\zeta ,t)$ and $I(\zeta ,t)$ results in the equations

which are 2*N* deterministic equations which are a function of the deterministic initial conditions $I0(0)$ and $VC0(0)$.

which results in 3*N* deterministic equations which are function of the initial conditions given in Eqs. (32*a*)–(32*c*).

*a*)–(27

*c*) can be rewritten as

#### 4.1.2 Numerical Results

*System model*$\mathbb{R}L\u21021$.

The three models $\mathbb{R}L\u21021,\u2009\mathbb{R}L\u21022$, and $\mathbb{R}L\u21023$ will be studied with $Vs=10\u2009\u2009V,\u2009Rl=2.7\u2009\u2009\Omega $, and $Ru=3.3\u2009\u2009\Omega $. The initial conditions $VL(t0)=3\u2009\u2009V\u2009and\u2009VC(t0)=4\u2009\u2009V$ in conjunction with the uncertain resistance make $I(t0)$ an uncertain variable. The coefficients of the PC expansion of $I(\zeta ,t0)$ are determined via Eq. (32*c*). 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*), $VC(t)$, and $VL(t)$ 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 $I(\zeta ,t)$ is nonzero at the initial time reflecting the uncertainty in the initial value of the current generated by the uncertainty in the resistance.

*System model*$\mathbb{R}L\u21022$.

In the attempt to reduce the index of the DAE, a simple method of substitution was used (as explained in Eqs. (28*a*) and (28*b*)) 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 $\mathbb{R}L\u21022$ is solved by using “$ode45$” function of matlab. Sixth-order PC expansion was used to represent the system states $VC(\zeta ,t)$ and $I(\zeta ,t)$.

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 $I(\zeta ,t)$ at the initial time was observed.

*System model*$\mathbb{R}L\u21023$.

In the attempt to reduce the index of the DAE (the method explained in Eqs. (29*a*)–(29*c*) 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. (37*a*)–(37*c*) 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., $I(t0,\zeta )$, 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 $\mathbb{R}L\u21021$ and $\mathbb{R}L\u21022$.

### 4.2 Simple Pendulum.

It can be observed that the DAEs are not of semi-explicit form. The algebraic constraint in Eq. (48*b*) needs to be differentiated two times in order to express it into a fully implicit ODE.

for all the instances of time *t*.

*a*)–(51

*d*) represented as

since *A* is a nonsingular matrix for every instance of time *t*.

#### 4.2.1 Treatment of the Initial Conditions.

*L*of the pendulum rod. The resulting PC expansion of the uncertain model parameters is

*a*) and (48

*b*) must satisfy the explicit algebraic constraint Eq. (48

*b*) and the hidden constraint Eq. (53) at every instance of time

*t*, special treatment of the initial conditions similar to the RLC problem is necessary. The constraint equations are

*a*) and (58

*b*) expressed in terms of the PC coefficients lead to

*L*is considered to be uniformly distributed ($L\u223cU(Ll,Lu)$)), and the Wigner-3

*j*symbols will be used here as well. Galerkin projection of Eqs. (59) and (60) onto the basis $\varphi k(\zeta )$ results in the equations

*δ*is the Kronecker delta. Given that the length of the pendulum rod

_{ij}*L*has uncertainty, the algebraic constraint associated with the length of the pendulum results in the coupled initial conditions. The stochastic ODEs representing the equations of motion can be rewritten as

*a*) and (63

*b*) are linear, Galerkin projection onto the basis function $\varphi k(\zeta )$ results in the equations

*c*) and (63

*d*) which are bilinear are reduced to

*L*is the length of the pendulum. Galerkin projection of the stochastic states onto the basis function yields

Unlike Eqs. (37*a*)–(37*c*) where Galerkin projection resulted in closed-form expressions, the same is not possible for Eqs. (69*a*) and (69*b*) which will be referred to as $\mathbb{P}2$. 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.

*M*samples are drawn from the distribution

*t*is used to solve the least squares problem

_{k}Note that the pseudo-inverse term $(PTP)\u22121\u2009\u2009PT$ needs to only be evaluated once. A similar equation is derived to determine the PC coefficients of the $\omega (t,\zeta )$.

*a*)–(51

*d*) can be rewritten as

The Gauss quadrature formula using *N* quadrature points ensures an exact evaluation of integrals of polynomials of order $2N\u22121$. Since the terms in Eqs. (74*c*) and (74*d*) have a highest order of $2N\u22122$, inner product on the *N* − 1 order basis functions would result in polynomials of order $3N\u22123$ which would only result in the Gauss quadrature formula approximating the inner product integrals.

#### 4.2.2 Numerical Results.

*a*)–(51

*d*)) and the other given by Eq. (47), which is a second-order ODE. The results of the simulation of the model $\mathbb{P}2$ are transformed to Cartesian space using the equations

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 × 10^{5} 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 × 10^{5} 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 MC method and surrogate PC model calculated at different instance of time | ||||
---|---|---|---|---|

State | ||||

Time | Z_{1} | Z_{2} | Z_{3} | Z_{4} |

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 | Z_{1} | Z_{2} | Z_{3} | Z_{4} |

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).