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) [1116]. 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 [1721] 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

It is a common practice to rearrange higher-order ODEs into a vector of first-order differential equations when developing numerical solutions. These are initial value problem of the form
(1)
where F and y are vectors. The vector of differential equations is called implicit if the state derivative y(t) cannot be isolated and represented in terms of the states y. If Eq. (1) can be rewritten in the explicit or canonical form
(2)
with the same state variables y(t), then the system of ODEs is referred to as explicit ODEs [4], and there are numerous algorithms for solving the initial value problem. If Eq. (1) can be rewritten as
(3a)
(3b)

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

The algebraic constraints along with the system of ODEs can appear in a more special form
(4a)
(4b)
which is referred to as semi-explicit DAEs. y(t),G(t,y(t),F(t,y(t) are all vector valued. In this work, we will focus on both the aforementioned forms of DAEs given by Eqs. (3a), (3b), (4a), and (4b). Equations (3a) and (3b) can be rewritten as
(5a)
(5b)

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 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 (4b) and (5b), hidden constraints need to be enforced. Every equation of the set of k1 equations derived after jth differentiation (j=1,2,k) of F(t,y(t)) is called a hidden constraint.

An important observation in case of DAEs is that the index does not only depend on the form of the DAE (Eqs. (5a) and (5b)) 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:
(6a)
(6b)
(6c)

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

  1. Setting y2(t)=1, Eq. (6c) yields y3(t)=t. After substituting the solution for y3(t) in Eq. (6a), the solution for y1(t) yields y1(t)=(t2/2)+y1(t0) (t0 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.

  2. Setting y2(t)=1, Eq. (6c) yields y1(t)=t. After differentiating Eq. (6c) 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.

Another way of explaining the dependence of an index of a system of DAEs to the solution of the algebraic variables is by replacing Eq. (6b) with the first derivative of Eq. (6b), which yields
(7)

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

The resulting system of DAE yields
(8)
Now, it can still be seen that the index of the system of the DAEs depends on the initial conditions y2(t)=0 and y2(t)=1. The solution corresponding to the initial condition y2(t)=0 results in an index one system of DAEs, and the second solution yields an index two system of DAEs.
For general DAE systems (1) as defined in Ref. [6], the index along a solution y(t) is the minimum number of time (say, p times) the existing DAEs need to be differentiated to solve for y(t) uniquely
(9)
Consider the following special case of semi-explicit DAEs, which is linear DAEs with time varying coefficients:
(10)
By differentiating the constraint Eq. (10) with respect to t, the result is
(11)

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.

Both Eqs. (3a) and (3b) and Eqs. (4a) and (4b) can be generalized to represent stochastic DAEs with probabilistic time-invariant model parameters characterizing the system. Equations (3a) and (3b) being the most general form, the polynomial chaos method to quantify uncertainty for that class of DAEs which are represented as
(12a)
(12b)

will be developed, where ζ is a random vector with probably density function (PDF) P(ζ). 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 WienerAskey 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.

If Φ=[ϕ0(ζ)ϕ1(ζ)ϕ2(ζ)] is the vector of the orthogonal polynomial functions corresponding to the PDF, P(ζ), the system states can be written as
(13)

For computational tractability, the PC expansion is truncated to a finite order which results in an approximation of the exact solution.

By taking Galerkin projection of Eqs. (12a) and (12b) on each of the orthogonal polynomial basis, Eq. (14) result
(14a)
(14b)

Hence, for N number of orthogonal bases, the number of initial equations (dim(G)+dim(g)) become (dim(G)+dim(g))×N number of equations (× stands for multiplication operation here), which are needed to be solved simultaneously in order to find dim(y(t,ζ))×N number of unknown functions of time.

The expressions for mean and variance can also be calculated using the method of intrusive PC. Since the set of orthogonal polynomials are created to be orthogonal with respect to a weighting function which is the probability density function of the random input, the mean is given as
(15)
Similarly, the second moment is given as
(16)
where
(17)

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.

For illustrative purposes, assume a stochastic model of the form
(18)
where ζU(1,1) is a uniformly distributed scalar, and xRn is the state vector of the system. A surrogate model using a PC expansion with Lagrange polynomial function will be used to develop a deterministic model. The Lagrange interpolation function
(19)
where ζi is a quadrature point, are orthogonal since they have an inner product
(20)

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 P(ζ). For example, for a uniformly distributed variable, ζi are the roots of the Legendre polynomial.

The Gauss quadrature formula states that there are positive weights w1,w2,,wn such that if f(x,ζ)Li(ζ) is a polynomial of degree at most 2n1, the equation
(21)

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.

A finite PC expansion of the states leads to
(22)
where XiRn, resulting in the model
(23)
Galerkin projection of Eq. (23) on the Lagrange polynomial Lm(ζ) yields the set of equations
(24a)
(24b)

for m=1,2,,N. If the terms on the left-hand and right-hand side of Eq. (18) are polynomials of degree less than or equal to 2n1, then Eq. (24b) represents Eq. (24a) exactly, else the approximation tends to the exact solution as n [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.

In this RLC circuit, Vs(t) is the control input voltage, and VR(t),VL(t),andVC(t) denote voltage across the resistor (R), the inductor (L), and the capacitor (C), respectively. Then, from Kirchoff's laws, the following circuit equations (DAEs) are derived (Fig. 1):
(25)
Fig. 1
Single loop RLC circuit diagram
Fig. 1
Single loop RLC circuit diagram
Close modal
Another way, the circuit equations can be written as follows:
(26a)
(26b)
(26c)
(26d)

In these DAEs, there are four system states (I(t),VL(t),VR(t),andVC(t)). In general, VsVs(t) ( denotes as equivalent to) is a function of time and is a known input.

Equation (26c) yields VRVR(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:
(27a)
(27b)
(27c)

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 L1 in this work.

The fully explicit form or canonical form of the system can be written in two ways. First, the equations can be written in terms of two states (for example, I(t) and VC(t)) as follows:
(28a)
(28b)

After solving the system of the ODEs, VR and VL can be computed from Eqs. (26c) and (27c), respectively. The system described by Eq. (28) will be referred to as L2 in this work.

In contrast to the ODEs given by Eqs. (28a) and (28b), the system model can be represented as an ODE by differentiating Eq. (27c) and substituting Eqs. (27a) and (27b) in the resulting differential equation, which is represented in terms of three states I(t),VC(t), and VL(t)) as
(29a)
(29b)
(29c)

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

For illustrative purposes, in this problem, the uncertain resistance is considered to be uniformly distributed (RU(Rl,Ru)). Since the uncertainty has a uniform distribution, the polynomial bases for the PC are Legendre polynomials. The uncertain resistance is represented as
(30)
and the PC expansion of all the uncertain variables is represented as
(31a)
(31b)
(31c)
(31d)

and ϕi(ζ) are Legendre polynomials, and ° denotes the dot product of two vectors.

Unlike in ODEs, the system states for the RLC circuit model (L1) are algebraically constrained at every instance of time by the equation
To illustrate the determination of the initial conditions for all the coefficients of the PC expansion of the system states, we assume that voltages VL(t0) and VC(t0) are independent initial conditions. This constrains the coefficients of the PC expansion of the state I(t0) at the initial time to be
(32a)
(32b)
(32c)
Galerkin projection is used to project Eq. (32c) onto the random space spanned by the polynomial basis ϕk(ζ). This requires triple product integrals with the weighting function P(ζ)=0.5. The identity
(33)
where the special case Wigner-3j symbol is
(34a)
(34b)
will be used to simplify Eq. (32c) to
(35)
where k=0,1,2,,(N1). To illustrate the consistent evolution of the uncertain states and their statistics, two equivalent models of the RLC circuit given by Eqs. (28a) and (28b) (L2) and Eqs. (29a)(29c) (L3) will be studied. Galerkin projection of the ODEs of model L2 after the polynomial expansion of VC(ζ,t) and I(ζ,t) results in the equations
(36a)
(36b)

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

The corresponding deterministic equations resulting from Galerkin projection of model L3 are
(37a)
(37b)
(37c)

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

Finally, a nonintrusive surrogate model using Lagrange interpolation polynomials as basis functions is used to characterize the statistics of the evolving states of the RLC circuit. Representing the states as
(38)
(39)
(40)
Equations (27a)(27c) can be rewritten as
(41)
(42)
(43)
Galerkin projection of Eqs. (41) and (42) on the Lagrange polynomial Lm(ζ) in conjunction with the Gauss quadrature formula leads to the equations
(44)
(45)
(46)

for m=1,2,,n. n simulations of Eqs. (44) and (45) in conjunction with initial conditions satisfying Eq. (46) permit us to determine the mean and variance of the evolving states of the RLC circuit.

4.1.2 Numerical Results

System modelL1.

The three models L1,L2, and L3 will be studied with Vs=10V,Rl=2.7Ω, and Ru=3.3Ω. The initial conditions VL(t0)=3VandVC(t0)=4V in conjunction with the uncertain resistance make I(t0) an uncertain variable. The coefficients of the PC expansion of I(ζ,t0) 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), 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(ζ,t) is nonzero at the initial time reflecting the uncertainty in the initial value of the current generated by the uncertainty in the resistance.

Fig. 2
Comparison of PC and Monte Carlo results from ℝLℂ1 model. (a) Estimate of state means and (b) estimate of state variances.
Fig. 2
Comparison of PC and Monte Carlo results from ℝLℂ1 model. (a) Estimate of state means and (b) estimate of state variances.
Close modal
System modelL2.

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 L2 is solved by using “ode45” function of matlab. Sixth-order PC expansion was used to represent the system states VC(ζ,t) and I(ζ,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(ζ,t) at the initial time was observed.

Fig. 3
Comparison of PC and Monte Carlo results from ℝLℂ2 model. (a) Estimate of state means and (b) estimate of state variances.
Fig. 3
Comparison of PC and Monte Carlo results from ℝLℂ2 model. (a) Estimate of state means and (b) estimate of state variances.
Close modal
System modelL3.

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., I(t0,ζ), 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 L1 and L2.

Fig. 4
Comparison of PC and Monte Carlo results from ℝLℂ3 model. (a) Estimate of state means and (b) estimate of state variances.
Fig. 4
Comparison of PC and Monte Carlo results from ℝLℂ3 model. (a) Estimate of state means and (b) estimate of state variances.
Close modal

4.2 Simple Pendulum.

The second benchmark problem considered is the single pendulum illustrated in Fig. 5 modeled in Cartesian space. The well recognized nonlinear governing equation of motion for a simple pendulum is
(47)
when represented in a polar coordinate system, but it is a DAEs' system when the governing equation of motion is represented in the Cartesian coordinate system
(48a)
(48b)
Fig. 5
Schematics of the simple pendulum
Fig. 5
Schematics of the simple pendulum
Close modal

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.

The resultant fully implicit ODEs are
(49a)
(49b)
The system can be represented in first-order form by defining the variables
(50)
and the implicit ODEs are rewritten as
(51a)
(51b)
(51c)
(51d)
The solution for z1,z2,z3, and z4 must also satisfy the following explicit algebraic constraint equation:
(52)
and the following “hidden constraint” equation resulting from the first derivative of the algebraic constraint Eq. (52):
(53)

for all the instances of time t.

Since the definition of fully implicit ODEs is that the system of ODEs cannot be expressed in standard canonical form of ODEs, Eqs. (51a)(51d) represented as
(54)
can be rewritten in the canonical form of the fully implicit system
(55)

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

4.2.1 Treatment of the Initial Conditions.

Similar to the benchmark RLC circuit problem, the single pendulum model is assumed to be characterized by uncertainty in one variable, the length L of the pendulum rod. The resulting PC expansion of the uncertain model parameters is
(56a)
(56b)
(56c)
(56d)
(56e)
where the augmented state space model includes as states the coefficients of the basis functions
(57)
Since the solution for the system states z1,z2,z3,andz4 from the system of DAE's Eqs. (48a) and (48b) must satisfy the explicit algebraic constraint Eq. (48b) 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
(58a)
(58b)
which implies that two initial conditions are dependent on the other two. Assuming that the initial conditions, z1(t0)(x(ζ,t0)x(t0)) and z3(t0)(x˙(ζ,t0)x˙(t0)), are independent initial conditions, Eqs. (58a) and (58b) expressed in terms of the PC coefficients lead to
(59)
(60)
A Galerkin projection is used to project Eqs. (59) and (60) onto the random space spanned by the polynomial basis. This requires triple product integrals with the weighting function P(ζ)=0.5 (in this problem, the uncertainty in L is considered to be uniformly distributed (LU(Ll,Lu))), and the Wigner-3j symbols will be used here as well. Galerkin projection of Eqs. (59) and (60) onto the basis ϕk(ζ) results in the equations
(61)
(62)
where δij is the Kronecker delta. Given that the length of the pendulum rod 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
(63a)
(63b)
(63c)
(63d)
Since Eqs. (63a) and (63b) are linear, Galerkin projection onto the basis function ϕk(ζ) results in the equations
(64)
(65)
and Eqs. (63c) and (63d) which are bilinear are reduced to
(66)
(67)
Equations (64)(67) which are implicit bilinear equations in conjunction with initial conditions given by the solution of Eqs. (61) and (62) are solved using the ode15i function of matlab. The resulting PC coefficients are used to characterize the mean and variance of the evolution of the uncertain states of the pendulum modeled as DAEs. The model will be referred to as 1. For the sake of completeness, a PC surrogate model has been developed for the pendulum modeled as ODEs with states (θ,θ˙) representing the angular displacement and velocity of the pendulum. The states are represented as
(68a)
(68b)
where L is the length of the pendulum. Galerkin projection of the stochastic states onto the basis function yields
(69a)
(69b)

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

Two nonintrusive approaches can be used to evaluate the mean and variance of the evolving states of the pendulum. The first is based on the least squares PC expansion. For the stochastic pendulum problem, M samples are drawn from the distribution
(70)
and Eq. (47) is simulated, and the time evolution of the states θ(tk,ζ),ω(tk,ζ) at discrete time instants tk is used to solve the least squares problem
(71)
whose solution is
(72)

Note that the pseudo-inverse term (PTP)1PT needs to only be evaluated once. A similar equation is derived to determine the PC coefficients of the ω(t,ζ).

The second nonintrusive PC expansion is based on the use of (N1)th order Lagrange interpolation polynomials as basis function. Representing the four states of the DAE representing the pendulum as
(73)
Equations (51a)(51d) can be rewritten as
(74a)
(74b)
(74c)
(74d)

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

4.2.2 Numerical Results.

To illustrate the estimation of time evolution of the mean and variance of the states of the pendulum represented in Cartesian coordinates which resulted in a set of DAE and in polar coordinates resulting in a set of ODEs, numerical simulations with initial conditions of x(t0)=z1(t0)=4m and x˙(t0)=z3(t0)=0m/s and an uncertain pendulum length with bounds of Ll=4.5m and Lu=5.5m are executed. Multiple order of PC expansions are studied to study convergence, and it is noted that when four or more basis functions are used, there is no noticeable change in estimates of the mean and variance of the states. Results presented in this section correspond to using six Legendre polynomials as basis functions, while 100,000 Monte Carlo simulations are used to determine the statistical mean and variance and will be used to benchmark the results of the PC surrogate models 1 and 2. It should be noted that two sets of Monte Carlo simulations are conducted; one using the implicit DAEs (Eqs. (51a)(51d)) and the other given by Eq. (47), which is a second-order ODE. The results of the simulation of the model 2 are transformed to Cartesian space using the equations
(75)
(76)
(77)
(78)

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.

Fig. 6
Comparison of means calculated by PC and MC
Fig. 6
Comparison of means calculated by PC and MC
Close modal

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.

Fig. 7
Comparison of variances calculated by PC and MC
Fig. 7
Comparison of variances calculated by PC and MC
Close modal

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.

Fig. 8
Evolution PDF of all states of DAEs (model ℙ1) using Monte Carlo simulations
Fig. 8
Evolution PDF of all states of DAEs (model ℙ1) using Monte Carlo simulations
Close modal
Fig. 9
Evolution PDF of all states of DAEs (model ℙ1) using polynomial chaos method
Fig. 9
Evolution PDF of all states of DAEs (model ℙ1) using polynomial chaos method
Close modal
Table 1

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
TimeZ1Z2Z3Z4
0.00 s9.9565 × 10−77.1807 × 10−500
0.45 s3.2802 × 10−75.0317 × 10−55.7881 × 10−79.0814 × 10−5
0.90 s1.0056 × 10−64.0612 × 10−63.3365 × 10−61.1582 × 10−4
1.35 s1.0185 × 10−63.4446 × 10−53.6109 × 10−69.5807 × 10−5
1.80 s3.2616 × 10−67.4939 × 10−51.583 × 10−51.4622 × 10−5
2.25 s7.062 × 10−65.1427 × 10−54.7912 × 10−61.0248 × 10−4
2.70 s6.9795 × 10−64.1967 × 10−65.9642 × 10−58.0802 × 10−5
3.15 s2.4678 × 10−52.2153 × 10−54.5262 × 10−53.2469 × 10−5
3.60 s1.8977 × 10−54.7335 × 10−59.7543 × 10−51.2413 × 10−4
4.05 s4.5417 × 10−53.0854 × 10−52.3585 × 10−57.7478 × 10−5
4.50 s3.1358 × 10−52.4963 × 10−51.0882 × 10−43.3921 × 10−4
Wasserstein distance between MC method and surrogate PC model calculated at different instance of time
    State
TimeZ1Z2Z3Z4
0.00 s9.9565 × 10−77.1807 × 10−500
0.45 s3.2802 × 10−75.0317 × 10−55.7881 × 10−79.0814 × 10−5
0.90 s1.0056 × 10−64.0612 × 10−63.3365 × 10−61.1582 × 10−4
1.35 s1.0185 × 10−63.4446 × 10−53.6109 × 10−69.5807 × 10−5
1.80 s3.2616 × 10−67.4939 × 10−51.583 × 10−51.4622 × 10−5
2.25 s7.062 × 10−65.1427 × 10−54.7912 × 10−61.0248 × 10−4
2.70 s6.9795 × 10−64.1967 × 10−65.9642 × 10−58.0802 × 10−5
3.15 s2.4678 × 10−52.2153 × 10−54.5262 × 10−53.2469 × 10−5
3.60 s1.8977 × 10−54.7335 × 10−59.7543 × 10−51.2413 × 10−4
4.05 s4.5417 × 10−53.0854 × 10−52.3585 × 10−57.7478 × 10−5
4.50 s3.1358 × 10−52.4963 × 10−51.0882 × 10−43.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).

References

1.
Musau
,
P.
,
Lopez
,
D. M.
,
Tran
,
H.-D.
, and
Johnson
,
T. T.
,
2018
, “
Linear Differential-Algebraic Equations (Benchmark Proposal)
,”
EPiC Ser. Comput.
,
54
, pp.
174
184
.10.29007/4gj7
2.
Campbell
,
S.
,
Ilchmann
,
A.
,
Mehrmann
,
V.
, and
Reis
,
T.
,
2019
,
Applications of Differential-Algebraic Equations: Examples and Benchmarks
,
Springer
, Cham, Switzerland.
3.
Lamour
,
R.
,
März
,
R.
, and
Weinmüller
,
E.
,
2015
, “
Boundary-Value Problems for Differential-Algebraic Equations: A Survey
,”
Surveys in Differential-Algebraic Equations III
,
Springer
, Cham, Switzerland, pp.
177
309
.
4.
Brenan
,
K. E.
,
Campbell
,
S. L.
, and
Petzold
,
L. R.
,
1995
,
Numerical Solution of Initial-Value Problems in Differential-Algebraic Equations
,
SIAM
, Philadelphia, PA.
5.
März
,
R.
,
2015
, “
Differential-Algebraic Equations From a Functional-Analytic Viewpoint: A Survey
,”
Surveys in Differential-Algebraic Equations II
(Differential-Algebraic Equations Forum),
A.
Ilchmann
and
T.
Reis
, eds., Springer, Cham, Switzerland, pp.
163
285
.
6.
Ascher
,
U. M.
, and
Petzold
,
L. R.
,
1998
,
Computer Methods for Ordinary Differential Equations and Differential-Algebraic Equations
, Vol.
61
,
SIAM
, Philadelphia, PA.
7.
Brown
,
P. N.
,
Hindmarsh
,
A. C.
, and
Petzold
,
L. R.
,
1994
, “
Using Krylov Methods in the Solution of Large-Scale Differential-Algebraic Systems
,”
SIAM J. Sci. Comput.
,
15
(
6
), pp.
1467
1488
.10.1137/0915088
8.
Brown
,
P. N.
,
Hindmarsh
,
A. C.
, and
Petzold
,
L. R.
,
1998
, “
Consistent Initial Condition Calculation for Differential-Algebraic Systems
,”
SIAM J. Sci. Comput.
,
19
(
5
), pp.
1495
1512
.10.1137/S1064827595289996
9.
Shampine
,
L. F.
,
2002
, “
Solving 0 = f (t, y(t), dy/dt(t)) in Matlab
,”
J. Numer. Math.
,
10
(
4
), pp.
291
310
.10.1515/JNMA.2002.291
10.
Kunkel
,
P.
, and
Mehrmann
,
V.
,
2006
,
Differential-Algebraic Equations: Analysis and Numerical Solution
, Vol.
2
,
European Mathematical Society
, Zurich, Switzerland.
11.
Ghanem
,
R.
, and
Spanos
,
P. D.
,
1990
, “
Polynomial Chaos in Stochastic Finite Elements
,”
ASME J. Appl. Mech.
,
57
(
1
), pp.
197
202
.10.1115/1.2888303
12.
Kim
,
K.-K. K.
,
Shen
,
D. E.
,
Nagy
,
Z. K.
, and
Braatz
,
R. D.
,
2013
, “
Wiener's Polynomial Chaos for the Analysis and Control of Nonlinear Dynamical Systems With Probabilistic Uncertainties [Historical Perspectives]
,”
IEEE Control Syst. Mag.
,
33
(
5)
, pp.
58
67
.10.1109/MCS.2013.2270410
13.
Eldred
,
M.
, and
Burkardt
,
J.
,
2009
, “
Comparison of Non-Intrusive Polynomial Chaos and Stochastic Collocation Methods for Uncertainty Quantification
,”
AIAA
Paper No. 2009–976.
14.
Wiener
,
N.
,
1938
, “
The Homogeneous Chaos
,”
Am. J. Math.
,
60
(
4
), pp.
897
936
.10.2307/2371268
15.
Le Maître
,
O.
, and
Knio
,
O. M.
,
2010
,
Spectral Methods for Uncertainty Quantification: With Applications to Computational Fluid Dynamics
,
Springer Science & Business Media
, Dordrecht, The Netherlands.
16.
Xiu
,
D.
,
2010
,
Numerical Methods for Stochastic Computations: A Spectral Method Approach
,
Princeton University Press
, Princeton, NJ.
17.
Pulch
,
R.
,
2009
, “
Polynomial Chaos for Multirate Partial Differential Algebraic Equations With Random Parameters
,”
Appl. Numer. Math.
,
59
(
10
), pp.
2610
2624
.10.1016/j.apnum.2009.05.015
18.
Pulch
,
R.
,
2008
, “
Polynomial Chaos for Analysing Periodic Processes of Differential Algebraic Equations With Random Parameters
,”
PAMmml: Proceedings in Applied Mathematics and Mechanics
, Vol.
8
,
Wiley Online Library
, Weinheim, Germany, pp.
10069
10072
.
19.
Pulch
,
R.
,
2011
, “
Polynomial Chaos for Linear Differential Algebraic Equations With Random Parameters
,”
Int. J. Uncertainty Quantif.
,
1
(
3
), pp.
223
240
.10.1615/Int.J.UncertaintyQuantification.v1.i3.30
20.
Pulch
,
R.
,
2013
, “
Polynomial Chaos for Semiexplicit Differential Algebraic Equations of Index 1
,”
Int. J. Uncertainty Quantif.
,
3
(
1
), pp.
1
23
.10.1615/Int.J.UncertaintyQuantification.2011003306
21.
Pulch
,
R.
,
2014
, “
Stochastic Collocation and Stochastic Galerkin Methods for Linear Differential Algebraic Equations
,”
J. Comput. Appl. Math.
,
262
, pp.
281
291
.10.1016/j.cam.2013.10.046
22.
Dai
,
L.
,
1989
,
Singular Control Systems
(Lecture Notes in Control and Information Sciences), Vol.
118
, Springer, Berlin/Heidelberg, Germany.
23.
Xiu
,
D.
, and
Karniadakis
,
G. E.
,
2002
, “
The Wiener–Askey Polynomial Chaos for Stochastic Differential Equations
,”
SIAM J. Sci. Comput.
,
24
(
2
), pp.
619
644
.10.1137/S1064827501387826
24.
Lefebvre
,
T.
,
2021
, “
On Moment Estimation From Polynomial Chaos Expansion Models
,”
IEEE Control Syst. Lett.
,
5
(
5
), pp.
1519
1524
.10.1109/LCSYS.2020.3040851
25.
Nandi
,
S.
, and
Singh
,
T.
,
2019
, “
Nonintrusive Global Sensitivity Analysis for Linear Systems With Process Noise
,”
ASME J. Comput. Nonlinear Dyn.
,
14
(
2
), p.
021003
.10.1115/1.4041622
26.
Andrews
,
G. E.
,
Askey
,
R.
, and
Roy
,
R.
,
1999
,
Special Functions
, Vol.
71
,
Cambridge University Press
, Cambridge, UK.