Доступ предоставлен для: Guest
International Journal for Uncertainty Quantification
Главный редактор: Habib N. Najm (open in a new tab)
Ассоциированный редакторs: Dongbin Xiu (open in a new tab) Tao Zhou (open in a new tab)
Редактор-основатель: Nicholas Zabaras (open in a new tab)

Выходит 6 номеров в год

ISSN Печать: 2152-5080

ISSN Онлайн: 2152-5099

The Impact Factor measures the average number of citations received in a particular year by papers published in the journal during the two preceding years. 2017 Journal Citation Reports (Clarivate Analytics, 2018) IF: 1.7 To calculate the five year Impact Factor, citations are counted in 2017 to the previous five years and divided by the source items published in the previous five years. 2017 Journal Citation Reports (Clarivate Analytics, 2018) 5-Year IF: 1.9 The Immediacy Index is the average number of times an article is cited in the year it is published. The journal Immediacy Index indicates how quickly articles in a journal are cited. Immediacy Index: 0.5 The Eigenfactor score, developed by Jevin West and Carl Bergstrom at the University of Washington, is a rating of the total importance of a scientific journal. Journals are rated according to the number of incoming citations, with citations from highly ranked journals weighted to make a larger contribution to the eigenfactor than those from poorly ranked journals. Eigenfactor: 0.0007 The Journal Citation Indicator (JCI) is a single measurement of the field-normalized citation impact of journals in the Web of Science Core Collection across disciplines. The key words here are that the metric is normalized and cross-disciplinary. JCI: 0.5 SJR: 0.584 SNIP: 0.676 CiteScore™:: 3 H-Index: 25

Indexed in

ITERATIVE METHODS FOR SCALABLE UNCERTAINTY QUANTIFICATION IN COMPLEX NETWORKS

Amit Surana

United Technologies Research Center, East Hartford, Connecticut 06118, USA

Tuhin Sahai

United Technologies Research Center, East Hartford, Connecticut 06118, USA

Andrzej Banaszuk

United Technologies Research Center, East Hartford, Connecticut 06118, USA

Abstract

In this paper we address the problem of uncertainty management for robust design, and verification of large dynamic networks whose performance is affected by an equally large number of uncertain parameters. Many such networks (e.g., power, thermal, and communication networks) are often composed of weakly interacting subnetworks.We propose intrusive and nonintrusive iterative schemes that exploit such weak interconnections to overcome the dimensionality curse associated with traditional uncertainty quantification methods (e.g., generalized polynomial chaos, probabilistic collocation) and accelerate uncertainty propagation in systems with a large number of uncertain parameters. This approach relies on integrating graph theoretic methods and waveform relaxation with generalized polynomial chaos, and probabilistic collocation, rendering these techniques scalable. We introduce an approximate Galerkin projection that based on the results of graph decomposition computes “strong” and “weak” influence of parameters on states. An appropriate order of expansion, in terms of the parameters, is then selected for the various states. We analyze convergence properties of this scheme and illustrate it in several examples.

KEYWORDS: collocation, polynomial chaos, graph decomposition, waveform relaxation


1. INTRODUCTION

The issue of management of uncertainty for robust system operation is of interest in a large family of complex networked systems. Examples include power, thermal, and communication networks which arise in several instances such as more electric aircrafts, integrated building systems, and sensor networks. Such systems typically involve a large number of heterogeneous, connected components, whose dynamics is affected by possibly an equally large number of uncertain parameters and disturbances.

Uncertainty quantification (UQ) methods provide means of calculating probability distribution of system outputs, given the probability distribution of input parameters. Outputs of interest could include, for example, latency in communication network, power quality and stability of power networks, and energy usage in thermal networks. The standard UQ methods such as Monte Carlo (MC) [1] exhibit poor convergence rates whereas others such as quasi Monte Carlo (QMC) [2, 3], generalized polynomial chaos (gPC) [4], and the associated probabilistic collocation method (PCM) [5] suffer from the curse of dimensionality (in parameter space), and become practically infeasible when applied to a network as a whole. Improving these techniques to alleviate the curse of dimensionality is an active area of current research (see [6] and references therein): notable methods include the sparse-grid collocation method [7, 8] and analysis of variance (ANOVA) decomposition [9] for sensitivity analysis and dimensional reduction of the uncertain parametric space. However, none of such extension exploits the underlying structure and dynamics of the networked systems. In fact, many networks of interest (e.g., power, thermal, and communication networks) are often composed of weakly interacting subsystems. As a result, it is plausible to simplify and accelerate the simulation, analysis, and uncertainty propagation in such systems by suitably decomposing them. For example, authors in [10, 11] used graph decomposition to facilitate stability and robustness analysis of large-scale interconnected dynamical systems. Mezic et al. [12] used graph decomposition in conjunction with Perron Frobenius operator theory to simplify the invariant measure computation and uncertainty quantification, for a particular class of networks. While these approaches exploit the underlying structure of the system, they do not take advantage of the weakly coupled dynamics of the subsystems.

In this paper, we propose an iterative UQ approach that exploits the weak interactions among subsystems in a networked system to overcome the dimensionality curse associated with traditional UQ methods. We refer to this approach as probabilistic waveform relaxation (PWR), and propose both intrusive and nonintrusive forms of PWR. PWR relies on integrating graph decomposition techniques and waveform relaxation scheme, with gPC and PCM. Graph decomposition to identify weakly interacting subsystems can be realized by spectral graph theoretic techniques [13, 14]. Waveform relaxation [15] (WR), a parallelizable iterative method exploits this decomposition and evolves each subsystem forward in time independently but coupled with the other subsystems through their solutions from the previous iteration. In the intrusive PWR, the subsystems obtained from decomposing the original system are used to impose a decomposition on system obtained by Galerkin projection based on the gPC expansion. Further, the weak interactions are used to discard terms which are expected to be insignificant in the gPC expansion, leading to what we call an approximate Galerkin projected (AGP) system. We then propose to apply WR relaxation on the decomposed AGP system to accelerate the UQ computation. In the nonintrusive form of PWR, rather than deriving the AGP system, one works directly with subsystems obtained from decomposing the original system. At each waveform relaxation iteration we propose to apply PCM at the subsystem level, and use gPC to propagate the uncertainty among the subsystems. Bevause UQ methods are applied to relatively simpler subsystems which typically involve a few parameters, this renders a scalable nonintrusive iterative approach to UQ.We prove convergence of the PWR approach under very general conditions. Note that spectral graph decomposition can be done completely in a distributed fashion using a recently developed wave equation-based clustering method [16]. Moreover, one can further exploit timescale separation in the system to accelerate WR using an adaptive form of WR [17]. PWR when combined with wave equation-based distributed clustering and adaptive WR can lead to a highly scalable and computationally efficient approach to UQ in complex networks. Note that to the best of the authors′ knowledge, it is the first attempt to extend polynomial chaos-based methods to parallel environments such that every step in the process, from clustering to the simulations, is completely decentralized.

This paper is organized in six sections. In Section 2 we give the mathematical preliminaries for setting up the UQ problem for networked dynamical systems, and present an overview of gPC and PCM techniques. In Section 3 we discuss graph decomposition and waveform relaxation methods, which form basic ingredients of PWR. Here we also describe adaptive WR and wave equation based distributed graph decomposition techniques. We introduce the intrusive and nonintrusive PWR in Section 4 through a simple example, and then describe these methods in a more general setting. We also prove convergence of PWR, and analyze the scalability of the method. In Section 5 we illustrate the intrusive and nonintrusive PWR in several examples. Finally, in Section 6 we summarize the main results of this paper, and present some future research directions.

2. UNCERTAINTY QUANTIFICATION IN NETWORKED SYSTEMS

Consider a nonlinear system described by a system of random differential equation

(1)

where f = (f1, f1, ..., fn) ∈ ℝn is a smooth vector field, x = (x1, x1, ..., xn) ∈ ℝn are state variables, ξi 2 Rpi is a vector of random variables affecting the ith system. Let ξ = (ξT1, ..., ξTn)T ∈ ℝp be the p = ∑ni=1 pi dimensional random vector of uncertain parameters affecting the complete system. The solution to the initial value problem x(t0) = x0 will be denoted by x(t, ξ), where for brevity we have suppressed the dependence of solution on initial time t0 and initial condition x0. We shall assume that the system (1) is Lipschitz

(2)

where the Lipschitz constant L(ξ) depends on the random parameter vector and || · || is a Euclidean norm. We will assume that supξ∈ℝp L(ξ) = L < ∞.

Let us also define a set of quantities

(3)

as observables or quantities of interests. The goal is to numerically establish the effect of input uncertainty of ξ on output observables z. Naturally, the solution for system (1) and the observables (3) are functions of the same set of random variables ξ, i.e.,

(4)

In what follows we will adopt a probabilistic framework and model ξ = (ξ1, ξ2, ..., ξp) as a p-variate random vector with independent components in the probability space (Ω,𝒜,𝒫), whose event space is Ω and is equipped with σ- algebra 𝒜 and probability measure 𝒫. Throughout this paper, we will assume that the parameters ∑ = {ξ1, ξ2, ..., ξp} are mutually independent of each other. Let ρi: Γi → ℝ+ be the probability density of the random variable ξi(ω), with Γi = ξi(Ω) ⊂ ℝ being its image. Then, the joint probability density of any random parameter subset Λ = {ξi1, ξi2 , ..., ξim} ⊂ ∑ is given by

(5)

with a support

(6)

where | · | denotes the cardinality of the set. Without loss of generality we will assume that Γi = [−1 1], i = 1, ..., p.

Remark 2.1. While throughout the paper we will work with ordinary differential equation (ODE) (1) with parametric uncertainty, the PWR framework developed in this paper can be naturally extended to deal with (1) system of differential algebraic equations (DAEs), and (2) time-varying uncertainty. Both these extensions are illustrated through examples in Section 5.

2.1 Uncertainty Quantification Methods

In this section, we describe two inter-related UQ approaches: generalized polynomial chaos and the probabilistic collocation method. The gPC is an intrusive approach which requires explicit access to system equations (1), while PCM is a related sampling-based nonintrusive [and hence treats the system (1) as a black box] way of implementing gPC.

2.1.1 Generalized Polynomial Chaos

In the finite dimensional random space Γ defined in (6), the gPC expansion seeks to approximate a random process via orthogonal polynomials of random variables. Let us define one-dimensional orthogonal polynomial space associated with each random variable ξk, k = 1, ..., p as

(7)

where {ψik)}dki=0 denotes the orthonormal polynomial basis from the so-called Wiener-Askey polynomial chaos [4]. The Askey scheme of polynomials contains various classes of orthogonal polynomials such that their associated weighting functions coincide with probability density function of the underlying random variables. The corresponding P-variate orthogonal polynomial space in Γ is defined as

(8)

where the tensor product is over all possible combinations of the multi-index d = (d1, d2, ..., d|∑|) ∈ ℕ|∑| in set ℙ,

(9)

and P = (P1, ..., P|∑|)T ∈ ℕ|∑| is a vector of integers which restricts the maximum order of expansion of the ith variable ξi to be Pi, and P = maxi Pi. Thus, W∑,P is the space of N-variate orthonormal polynomials of total degree at most P with additional constraints on individual degrees of polynomials, and its basis functions Ψ∑,Pi (ξ) satisfy

(10)

for all 1 ≤ i, jN = dim(W∑,P). Note that in standard gPC all expansion orders are taken to be identical, i.e., P1,= P2... = P|∑| = P, so that dim(W∑,p) = (P + |∑|)!/P!|∑|!. We will however take advantage of an adaptive expansion, a notion which will be fully developed in Section 4.3.

The major advantage of applying the gPC is that a random differential equation can be transformed into a system of deterministic equations. A typical approach is to employ a stochastic Galerkin projection, in which all the state variables are expanded in polynomial chaos basis with corresponding modal coefficients [aik(t)], as

(11)

where the sum has been truncated to a finite order. Substituting these expansions in Eq. (1), and using the orthogonality property of polynomial chaos (10), we obtain for k = 1, ..., n, j = 1, ..., N,

(12)

a system of deterministic ODEs describing the evolution of the modal coefficients, with initial conditions

(13)

and a = (a11, ..., aN∑1, ..., a1n, ..., aNn)T,

(14)

with x∑,P(t, ξ) = [x∑,P1(t, ξ), ..., x∑,Pn(t, ξ)]. This system can be solved with any numerical method dealing with initial-value problems, e.g., the Runge-Kutta method. Similarly, the observable can be expanded in the gPC basis as

(15)

where

(16)

with k = 1, ..., d. Hence, once the solution to the system (12) has been obtained, the coefficients bjk can be approximated as

(17)

Such Galerkin procedures have been used extensively in the literature. In many instances Galerkin projection may not be possible due to the unavailability of direct access to the system of equations (1). In many other instances such intrusive methods are not feasible even in cases when the system equations are available, because of the cost of deriving and implementing a Galerkin system within available computational tools. To circumvent this difficulty, the probabilistic collocation method has been developed.

2.1.2 Probabilistic Collocation Method

PCM is a nonintrusive approach to solving stochastic random processes with the gPC. Instead of projecting each state variable onto the polynomial chaos basis, the collocation approach evaluates the integrals of form (16) by evaluating the integrand at the roots of the appropriate basis polynomials [5]. Given a one-dimensional (1D) probability density function ρjj), the PCM based on Gauss quadrature rule approximates an integral of a function g with respect to density ρjj) as follows:

(18)

where, rljkClj is the set of Gauss collocation points with associated weights wljk, lj is the accuracy level of quadrature formula, and mlj is the number of quadrature points corresponding to this accuracy level. Building on the 1D quadrature formula, the full grid PCM leads to following cubature rule:

(19)

where rlj = (rl1j1, ..., rlpjp), with l = (l1, ..., lp), and j = (j1, ..., jp) and wlj = wl1j1 ... wlnjn. To compute 𝒥(l, p) we need to evaluate the function on the full collocation grid C(l, p) which is given by the tensor product of 1D grids

(20)

with the total number of collocation points being Q = ∏pj=1 mlj . In this framework, therefore, for any t the approximations to the model coefficients ajk(t) [see Eq. (11)] and bjk(t) [see Eq. (15)] can be obtained as

(21)

with a similar expression for bjk(t). Note that to compute the summations arising in (21), the solution x(t, rlj) of the system (1) is required for each collocation point rlj in the full collocation grid C(l, p). Thus, the simplicity of collocation framework only requires repeated runs of deterministic solvers, without explicitly requiring the projection step in gPC.

If we choose the same order of collocation points in each dimension, i.e., ml1 = ml2, ... = mlpl, the total number of points is Q = lp, and the computational cost increases rather steeply with the number of uncertain parameters p. Hence, for large systems (n >> 1) with a large number of uncertain parameters (p >> 1), PCM becomes computationally intensive. As discussed in the Introduction, alleviating this curse of dimensionality is an active area of current research [6]. In this paper we propose a new uncertainty quantification approach which exploits the underlying network structure and dynamics to overcome the dimensionality curse associated with PCM. The key methodologies for accomplishing this are the graph decomposition and waveform relaxation, which are discussed in subsequent sections.

3. GRAPH DECOMPOSITION AND WAVEFORM RELAXATION

3.1 Waveform Relaxation

In this section we describe the basic mathematical concept of the waveform relaxation (WR) method for iteratively solving the system of differential equations of the form (1). For purposes of discussion here, we fix the parameter values ξ in the system (1) to fixed mean values. The general structure of a WR algorithm for analyzing system (1) over a given time interval [0, T] consists of two major steps: the assignment partitioning process and the relaxation process [15, 18].

Assignment-partitioning process: Let 𝒩 = {1, ..., n} be the set of state indices, and 𝒞i, i = 1, ..., m be a partition of 𝒩 such that

(22)

We shall represent by 𝒟: 𝒩𝓜 ≡ {1, 2, ..., m} a map which assigns the state index to its partition index, i.e., 𝒟(i) = j, where j is such that i𝒞j. Without loss of generality, we can rewrite Eq. (1) after the assignmentpartitioning process as

(23)

where, for each i = 1, ..., i,

(24)
(25)

with initial condition

(26)

and

(27)

are the subvectors assigned to the ith partitioned subsystem, such that ji k𝒞i; k = 1, ..., Mi = |𝒞i|, and

(28)

is a decoupling vector, with jk𝓜i and k = 1, ..., Ni = |𝓜i|. Here, 𝓜i is the set of indices of the partitions (or subsystems) with which the ith partition (or subsystem) interacts, and is given by

(29)

where 𝒩c i = {j𝒩: ∂Fi/∂xj = 0} and 𝔍(·) denotes the image of the map 𝒟.

Relaxation process: The relaxation process is an iterative procedure, with the following steps:

Step 1: (Initialization of the relaxation process) Set I = 1 and guess an initial waveform {y0i(t): t ∈ [0 T]} such that y0i(0) = y0i, i = 1, ..., m.

Step 2: (Analyzing the decomposed system at the Ith WR iteration) For each i = 1, ..., m, set

(30)

and solve the subsystem

(31)

over the interval [0, Ts] with initial condition yI i(0) = y0i, to obtain {yI (t): t ∈ [0, Ts]}.

Step 3: Set I = I + 1 and go to step 2 until satisfactory convergence is achieved.

Note that unlike traditional parallel DAE solvers, WR is agnostic to the numerical scheme [15, 19]. One can partition the original system based on timescales and then use appropriate numerical schemes (with different time steps) on various subsystems. This is particularly attractive for stiff dynamic systems [15]. Moreover, competing DAE solvers require communication at every function evaluation; this can be expensive for implicit methods since function evaluations are required for every Newton iteration [19]. WR requires communication only at the end of the simulation for the entire time window. The overall data communicated in WR also tend to be smaller since one is not forced to exchange information related to every function evaluation [20]. Since sending larger messages less frequently is more efficient from a communication point of view [21], the communication overhead of WR is significantly lower [19]. For a comparison of WR to competing approaches, see [15, 19, 20]. WR has been used in applications related to neural modeling [19], communication networks [17], and power grids [22]. Note that co-simulation [23] strategies for simulating heterogeneous systems, and utilizing different solvers for various subsystems, also frequently employ WR.

The general conditions for convergence of WR for a system of differential algebraic equations (DAEs) can be found in [17, 18]. Here, we recall a result from [17] specializing it for a system of differential equations.

Proposition 3.1. Convergence of WR for ODEs (see [18] for proof): Given that the system (1) is Lipschitz (condition 2), then for any initial piecewise continuous waveform {y0i(t): t ∈ [0, Ts]} such that y0i(0) = y0i [see definition (26)], i = 1, ..., m, the WR algorithm converges to the solution of (1) with initial condition x0.

A more intuitive analysis of error at each waveform iteration is described in [17]. Let y be the exact solution of the differential equation (23) and define EI to be the error of the Ith iterate, that is

(32)

As shown in [17], the error |EI| on the interval [0, T] is bounded as follows:

(33)

with C = eμT. Here C and η are related to the Lipschitz constants of the waveform relaxation operator [17]. Note here that the I! in the denominator dominates the error. Thus, with enough iterations one can make the error fall below any desired threshold. It is also evident from Eq. (33) that the error of standard waveform relaxation crucially depends on T. The longer the time interval, the greater is the number of iterations needed to bound the error below a desired tolerance. Based on this observation, a novel “adaptive” version of waveform relaxation has been developed in [17], which we refer to as adaptive waveform relaxation (AWR). The idea here is to perform waveform relaxation in “windows” of time that are picked so as to reduce I in Eq. (33). Specifically, one can pick small time intervals for computation when the solution to (1) changes significantly [implying E0(t) is large] and pick large intervals when the solution changes little [consequently E0(t) is small]. The solution from one time interval is extrapolated to the next using a standard extrapolation formula [24] and the initial error is estimated using

(34)

where

(35)

Here t0, t1, ..., tl are points through which one passes the extrapolating polynomial [17]. Note that φ(l) = dlφ/dtl is the lth derivative of the waveform relaxation operator φ with respect to t (see [17, 24]). We point out to the reader that upper bounds on the window size can be placed using the memory restrictions at the various computer processors. Also, at a particular processor, one does not need to store the waveforms of all the state variables. At a processor, one only needs to store the waveforms of state variables that influence the ones that are being computed locally.

The algorithm to compute the length of the time windows is as follows:

Adaptive waveform relaxation: To compute the time interval for ΔTI+1, execute the following steps:

  1. Set ΔTI+1 = 2ΔTI and δ = (1/20)ΔTI .
  2. Evaluate I+1,0TI + ΔTI+1) using Eq. (34) to estimate ||EI+1,0|| and compute ||ÊI+1,r|| with the aid of the following equation:
(36)
  1. If ||ÊI+1,r|| > ε and ΔTI+1 > (1/2)ΔTI , set ΔTI+1 = ΔTI+1 − δ and repeat step 2.

We define the minimal window length to be ΔTI+1 = (1/50)T. The above procedure gives a sequence of time intervals [0, T1], [T1, T2], ... , [Tν−1, Tν], where Tν = T, on which WR is performed (as described in relaxation process) with an initial “guess” waveform provided by an extrapolation of the solution on the previous interval [17]. AWR is found to accelerate simulations by orders of magnitude over traditional WR [17]. In this work, we propose to use AWR for simulating the set of differential equations that arise from intrusive polynomial chaos. As mentioned before, the curse of dimensionality gives rise to a combinatorial number of equations [4] making AWR particularly attractive.

While the convergence of WR or AWR is guaranteed irrespective of how the system is decomposed in the assignment-partitioning step, the rate of convergence depends on the decomposition [17]. For a given nonlinear system, determining a decomposition that leads to an optimal rate of AWR convergence is an NP-complete problem [17]. Ideally, to minimize the number of iterations required for convergence, one would like to place strongly interacting equations/variables on a single processor, with weak interactions between the variables or equations on different processors. In Eq. (33), η is a measure of the “strength” of connection between clusters. The smaller η is, the weaker are the interactions and consequently, fewer are the required iterations for convergence. In fact, once the system has been decomposed, by bounding η one can estimate the computational effort required by AWR, thus aiding in the decision of using this approach over sampling based methods.

We show in [17] that spectral clustering [13] along with horizontal vertical decomposition [12] is a good heuristic for decomposing systems for fast convergence in WR and AWR. For a comparison of spectral clustering with other approaches see [17]. For this task, we now discuss a novel decentralized spectral clustering approach [16] that when coupled with AWR [17] provides a powerful tool for simulating large dynamic systems, making every step of the UQ approach scalable.

3.2 Graph Decomposition

The problem of partitioning the system of equations (1) into subsystems based on how they interact or are coupled to each other can be formulated as a graph decomposition problem. Given the set of states x1, ..., xn and some notion of dependence wij ≥ 0, i = 1, ..., n, j = 1, ..., n between pairs of states, an undirected graph G = (V,E) can be constructed. The vertex set V = {1, ..., n} in this graph represent the states xi, and the edge set is EV × V, where a weight wij ≥ 0 is associated with each edge (i, j) ∈ E, and W = [wij] is the n × n weighted adjacency matrix of 𝒢. In order to quantify coupling strength wij between nodes or states, we propose to use

(37)

where J = {(1/Ts) ∫t0+Tst0 Jij [x(t, ξm), ξm, t]dt}, is the time average of the Jacobian,

(38)

computed along the solution x(t, ξ) of the system (1) for nominal values of parameters ξm. Use of the system Jacobian for horizontal vertical graph decomposition can also be found in [12].

We will now discuss spectral clustering (see [13], a popular graph decomposition/clustering approach that allows one to partition a undirected graph given its adjacency matrix W). In this method, first a (normalized) graph Laplacian is constructed as follows [14, 25, 26]:

(39)

or equivalently as L = ID−1W, where D is the diagonal matrix with the row sums of W. The clustering assignment/ decomposition is then obtained by computing the eigenvectors/eigenvalues of L. In particular, one uses the signs of the components of the second and higher eigenvectors to partition the nodes in the graph into clusters [13]. Traditionally, one can use standard matrix algorithms for eigenvector/eigenvalue computation [27]. However, as the size of the dynamic system or network (and thus corresponding adjacency matrix) increases, the execution of these standard algorithms becomes infeasible on monolithic computing devices. To address this issue, distributed eigenvalue/ eigenvector computation methods have been developed; see for example [28].

In [16], a wave equation-based distributed algorithm to partition large graphs has been developed which computes the partitions without constructing the entire adjacency matrix W of the graph [13]. In this method one “hears” clusters in the graph by computing the frequencies [using fast Fourier transform (FFT) locally at each node] at which the graph “resonates.” In particular, one can show that these “resonant frequencies” are related to the eigenvalues of the graph Laplacian L (39) and the coefficients of FFT expansion are the components of the eigenvectors. In fact, the algorithm is provably equivalent to the standard spectral clustering [13]; for details see [16].

The steps of the wave equation-based clustering algorithm are as follows. One starts by writing the local update equation at node i based on the discretized wave equation,

(40)

where ui(t) is the value of u at node i at time t and Lij are the local entries of the graph Laplacian. At each node i, ui(0) = ui(−1) is set to a random number on the interval [0; 1]. One then updates the value of ui using Eq. (40) until t = Tmax (for a discussion on how to pick Tmax see [16]). Note that ui(t) is a scalar quantity and one only needs nearest neighbor information in Eq. (40) to compute it. One then performs a local FFT on [ui(1), ..., ui(Tmax)] and then assigns the coefficients of the peaks of the FFT to vij. Here the frequency of the jth peak is related to λj , the jth eigenvalue of the graph Laplacian L, and vij is the ith component of the jth eigenvector.

In [16], it has been shown that for wave speed c < √2 in Eq. (40), the above wave equation-based algorithm is stable and converges. Moreover, the algorithm converges in O(√τ) steps, where τ is the mixing time of the Markov chain [16] associated with the graph G. The competing state-of-the-art algorithm [28] converges in O[τ(log(n)2)]. For large graphs or datasets, O(√τ) convergence is shown to provide orders of magnitude improvement over algorithms that converge in O[τ(log(n)2)]. For detailed analysis and derivations related to the algorithm, we refer the reader to [16].

4. ITERATIVE UNCERTAINTY QUANTIFICATION APPROACH

In this section we discuss how gPC and PCM can be integrated with the WR scheme, extending it to a probabilistic setting. As mentioned earlier, we refer to this iterative UQ approach as PWR. Figure 1 shows the schematic of PWR framework. In the intrusive PWR, the subsystems obtained from decomposing the original system are used to impose a decomposition on the system obtained by Galerkin projection based on the gPC expansion. Further, the weak interactions are used to discard terms which are expected to be insignificant in the gPC expansion, leading to what we call an approximate Galerkin projected (AGP) system. We then propose to apply standard or adaptive WR on the decomposed AGP system to accelerate the UQ computation. In the nonintrusive form of PWR, rather than deriving the AGP system, one works directly with subsystems obtained from decomposing the original system. At each waveform relaxation iteration we propose to apply PCM at subsystem level and use gPC to propagate the uncertainty among the subsystems. Note that unlike intrusive PWR (where deterministic decoupling vectors or deterministic waveforms are exchanged), in nonintrusive PWR stochastic decoupling vector or probabilistic waveforms represented in gPC basis are exchanged between subsystems at each iteration.

FIG. 1: Schematic of intrusive (left) and nonintrusive (right) PWR.

We first describe the key technical ideas behind intrusive and nonintrusive PWR though an illustration on a simple example in Section 4.1. These notions are fully generalized later in Sections 4.2–4.4. We also prove the convergence of the PWR approach (in Section 4.5), and in Section 4.6 discuss the computational gain it offers over standard application of gPC and PCM.

4.1 Main Ideas of PWR

We illustrate the proposed PWR framework through an example of parametric uncertainty in a simple system (1). Consider the following coupled oscillator system:

(41)

Here, ωi is the uncertain angular frequency of the ith (i = 1, 2, 3) oscillator. We assume that parameters ωi are mutually independent, each having probability density ρi = ρii). The coupling matrix K

(42)

is assumed deterministic with Kij = 𝒪(𝛜); 𝛜 << 1, so that the three oscillators in (1) weakly interact with each other, i.e., the subsystem 2 weakly affects subsystem 1 and 3, and vice versa.

4.1.1 Approximate Galerkin Projection for the Simple Example

In standard gPC, states xi are expanded in a polynomial chaos basis as

(43)

where Ψ∑,Pj η W∑,P, the P variate polynomial chaos space formed over ∑ = {ω123}, and P = (P1, P2, P3) determines the expansion order (see Section 2.1.1 for details). Note that in this expansion (43), the system states are expanded in terms of all the random variables ω affecting the entire system. From the structure of system (1) it is clear that the first subsystem is directly affected by the parameter ω1 and indirectly by parameter ω2 through the the state x2. We neglect the second-order effect of ω3 on x1. A similar statement holds true for subsystem 3, while subsystem 2 will be weakly influenced by ω1 and ω3 through states x1 and x3, respectively. This structure can be used to simplify the Galerkin projection as follows. For x1 we consider the gPC expansion over ∑1 = Λ1 ∪ Λc1,

(44)

where

(45)

and P1 = (P11, P12). Note that since ω2 weakly affects x1, the order of expansion P12 can be chosen to be smaller compared to P11. Similarly, one can consider simplification for x3,P3 3(t13). For x2 following similar steps, we define

(46)

and P2 = (P21, P22, P23). By similar argument, one will choose P21, P23 much smaller than P22. We also introduce the following two projections associated with the state x2:

(47)

where i = 1, 3 and 〈·, ·〉 is the appropriate inner product on W∑,P (see Section 4.3 for details).With these expansions, and using standard Galerkin projection we obtain the following system of deterministic equations:

(48)

with appropriate initial conditions, where

and a = (a1, a2, a3)T, with ai = (ai1, ..., ai1iNi) and F = (F1, F2, F3)T with Fi = (Fi1, ..., F1NΛi). We will refer to (48) as an AGP system. The notion of AGP in a more general setting is described in Section 4.3.

4.1.2 Intrusive PWR Illustrated on the Simple Example

In intrusive PWR, after performing the AGP explicitly, the system (48) is decomposed as

(49)

where d1 = 𝒫2,1(a2), d2 = (a1, a3) and d3 = 𝒫2,1(a2) are the decoupling vectors [here we overloaded notation for 𝒫2,i(a2) to imply the coefficients in expansion (47)]. Note that the decomposition of system (48) is based on the decomposition of the original system (41). Adaptive or standard WR, described in Section 3.1, can then be applied to solve the decomposed system (49) iteratively. Since, the the system (49) is deterministic, deterministic waveforms or deterministic decoupling vectors di, i = 1, 2, 3 are exchanged in each WR iteration (see Fig. 1 for an illustration).

4.1.3 Nonintrusive PWR Illustrated on the Simple Example

In the nonintrusive form of PWR, rather than deriving the AGP, one works directly with subsystems obtained from decomposing the original system. The main idea here is to apply PCM at subsystem level at each PWR iteration, use gPC to represent the probabilistic waveforms, and iterate among subsystems using these waveforms. Recall that in standard PCM approach (Section 2.1.2), the coefficients aim(t) are obtained by calculating the integral

(50)

The integral above is typically calculated by using a quadrature formula and repeatedly solving the ith subsystem over an appropriate collocation grid 𝒞i(∑i) = 𝒞ii) × 𝒞ic i), where, 𝒞ii) is the collocation grid corresponding to parameters Λi (and let ls be the number of grid points for each random parameter in Λi), 𝒞ic i) is the collocation grid corresponding to parameters Λc i (and let lc be the number of grid points for each random parameter in Λc i). Since the behavior of the ith subsystem is weakly affected by the parameters λc i , we can take a sparser grid in Λc i dimension, i.e., lc < ls, as we took lower-order expansion for these random variables in Section 4.1.1. Below we outline key steps in nonintrusive PWR:

Step 1: (Initialization of the relaxation process with no coupling effect): Set I = 1; guess an initial waveform x0i (t) consistent with

initial condition. Set d1 1 = x02, d12 = (x01, x03), d13 = x02, and solve

(51)

with an initial condition x1i (0)=x0i (0) on a collocation grid 𝒞ii). Determine the gPC expansion xΛi,Pi,1i (t, ·) by computing the expansion

coefficients from the quadrature formula (50).

Step 2: (Initialization of the relaxation process, incorporating first level of coupling effect): Set I = 2 and let d21 = xΛ2,P2,12, d22 = (xΛ1,P1,11,

xΛ3,P3,13), d23 = xΛ2,P2,12 be the stochastic decoupling vectors. Solve

(52)

over a collocation grid 𝒞i(∑i) to obtain xi,Pi,2i (t, ·). From now on we shall denote the solution of the ith subsystem at Ith iteration by

xi,Pi,I i.

Step 3: (Analyzing the decomposed system at the Ith iteration): Set the decoupling vectors, dI1 = P2,1(x2,P2,I−1 2), dI2 = (x1,I−11,

x3,P3,I−13), dI3 = P2,3(x2,P2,I−12), and solve

(53)

over a collocation grid Ci(∑i) and obtain the expansion xi,Pi,Ii (t, ·).

Step 4: (Iteration): Set I = I + 1 and go to step 5 until satisfactory convergence has been achieved.

Note that in the above nonintrusive PWR, the decoupling vectors are stochastic and so at each iteration probabilistic waveforms are exchanged between subsystems (see Fig. 1 for an illustration). We next generalize the intrusive and nonintrusive PWR introduced above in the forthcoming sections.

4.2 Decomposition of Galerkin Projected System

We begin by revisiting the complete Galerkin system (12). To apply WR, recall that the first step is the assignmentpartitioning (see Section 3.1). There are two possible approaches for partitioning the complete Galerkin system. One can first split the original dynamic system (1), and then use this decomposition to partition the complete Galerkin projection (12) by assigning the model coefficients in (11) for each state to the cluster to which it is assigned while decomposing system (1). As previously explained in Section 3.2, the partitioning is performed by representing the dynamic system (1) as a graph with the symmetrized time-averaged Jacobian (37) as the weighted adjacency matrix. One can then apply the wave equation-based decentralized clustering algorithm outlined in Section 3.2.

Alternatively, one can perform this decomposition directly on the complete Galerkin projection (12). Let the symmetrized time-averaged Jacobian for the resulting system (12) be

(54)

where = {(1/Ts) ∫ t0+Tst0 Jij [a(t), t]dt} is the time average of the Jacobian,

(55)

computed along the solution a(t) of the system (12). This gives

(56)

Taylor expanding fk[x(ξ, t), ξk, t] locally gives

(57)

Thus, one expects to get similar results by performing clustering on the original system [in (1)] to that obtained based on complete Galerkin system (12). Because the dimensionality of system (1) is much lower than that of system (12), the first decomposition is less computationally challenging than the latter. In this work, we use the original system to determine the decomposition and use that to impose the partition of the Galerkin projection. Given the decomposition of system (12), one can use the WR or its adaptive form to simulate the system in a parallel fashion. However, before doing this one can further exploit the weak interaction between subsystems to reduce the dimensionality of the complete Galerkin system, as described in Section 4.1.1. We next describe this approximate Galerkin projection in a more general setting.

4.3 Approximate Galerkin Projection and Intrusive Probabilistic Waveform Relaxation

Recall that in the gPC expansion (11), all the system states are expanded in terms of random variables affecting the entire system. However, the ith subsystem in the decomposition [see Eq. (23)] is directly affected by the parameters Λi [see definition (27)] and indirectly by other parameters through the decoupling vector [see definition (28)]. We shall denote by

(58)

the set of parameters which indirectly affect the ith subsystem through the immediate neighbor interactions, and denote by

(59)

the set of parameters that directly and indirectly (through nearest neighbor interaction) affect the ith subsystem. Under the hypothesis that the ith subsystem is dynamically weakly coupled with its nearest neighbors, uncertainty in parameters Λic will weakly influence the states in the ith subsystem through the decoupling vector, while the uncertainty in parameters ∑\∑i can be neglected. To capture this effect, consider a P-variate space (analogous to the P-variate space introduced in Section 2.1.1)

(60)

formed over any random parameter subset Λ = {ξi1, ξi2, ..., ξin} ⊂ ∑. We denote the basis elements of Wλ,P by ΨΛ,Pi, i = 1, ..., NΛ, where NΛ = dim(WΛ,P). Note that for any Λ1 ⊂ Λ2 ⊂ ∑,

(61)

where W = {0} is the P-variate space corresponding to the empty set. Also, recall that P1, P2 are vectors which control the expansion order in gPC expansion. The inner product on W∑,P induces an inner product on WΛ,P as follows

(62)

for any X1(ξ), X2(ξ) η WΛ,P. Using this inner product, we introduce a projection operator

(63)

such that for any X(ξ) η WΛ2,P2

(64)

With respect to the given decomposition 𝒟 imposed on the system (see Section 3.1), we define a projection operator 𝒫i,j indexed by subsystem i and state xj

Remark 4.1. For any subsystem i, since the parameters Λci weakly affect it, we can adaptively select component of vector Pi = (Pi,j1, ..., Pi,j|∑i|) so that a lower-order expansion is used in components corresponding to random variables in Λci.

For any subsystem i and a vector x∑,P (ξ, t) = [x∑,Pk1t), ..., x∑,Pknt)] [where x∑,Pit) expansion (11)], we associate a vector valued projection operator as follows:

(65)

In terms of these operators, for any state xk an approximate Galerkin projected equation is defined as

(66)

where i = 𝒟(k) is the index of the subsystem to which the state k belongs. More precisely, the above system can be expressed as

(67)

where a = (a𝒟(1) 11, ..., a𝒟(1) N𝒟(1),1, ..., a𝒟(n) 1n, ..., aN𝒟(n),n)T

(68)

and j = 1, ..., Ni , k = 1, ..., n with i = 𝒟(k). Let F = (F 𝒟(1)11, ..., F 𝒟(1) N𝒟(1),1, ..., F 𝒟(n) 1n , ..., F 𝒟(n) N𝒟(n),n)T; then, the system (67) can be compactly written as

(69)

with appropriate initial condition [see expression 13] and will be referred to as the approximate Galerkin projection (AGP). Using this generalization of AGP system, it is straightforward to generalize the intrusive PWR introduced in Section 4.1.2.

4.3.1 Intrusive PWR Algorithm

In the intrusive PWR, one applies the WR to the AGP system. Based on the discussion in Section 4.2, the decomposition 𝒟 of the original system (1) is used to impose a decomposition on the system (69), leading to

(70)

for i = 1, ..., m, where

(71)
(72)

ki𝒞i, and di = (aTji, ..., aTjNi)T is the decoupling vector (recall notation from Section 3.1). One then follows the procedure for waveform relaxation or its adaptive version, as described in Section 3.1. Adaptive WR can lead to a significant increase in convergence of WR as demonstrated in [17], and would be illustrated later in the Section 5.

As discussed in Section 2.1.1, the projection step (66) can be very tedious and in some cases not possible. Hence, applying waveform relaxation directly to the system (70) may not be practical in many instances. In the next section, we describe an alternative nonintrusive approach using probabilistic collocation, which does not require the projection step (66) explicitly.

4.4 Nonintrusive Probabilistic Waveform Relaxation

In terms of the projection operator (65), we can rewrite each subsystem in (23) as

(73)

where di(t, ·) is the stochastic decoupling vector or probabilistic waveform,

(74)

where we have explicitly indicated the dependence on the parameters [see definition (28)]. Here, for any i = 1, ..., m, yi,Pi = (xi,Piji1, ..., xi,PijiMi)T, with

(75)

By definition the coefficients amjik (t) in the above expansion satisfy the system (67). These coefficients can be obtained by using the quadrature formula (21), by repeatedly solving the system (73) over an appropriate collocation grid C(l, ni),

(76)

where l = (o, m), 𝒞(o, ni) = C1o1 × ... × C1oni is the collocation grid corresponding to parameters λi, with ni = |λi|, o = (o1, ..., oni), and 𝒞(m, nci) = C1 m1 × ... × C1mnci is the collocation grid corresponding to parameters λci, with nci = |λci| and m = (m1, ..., mnci). For simplicity we take o1 = ... = oni = ls and m1 = ... = mnci = lc for i = 1, ..., m. Since, the behavior of the ith subsystem is weakly affected by the parameters λci through the decoupling vector, then consistent with remark (4.1) we can take

(77)

leading to an adaptive collocation grid for each subsystem.With this, we are ready to generalize the nonintrusive PWR approach introduced in Section 4.1.3.

4.4.1 Non-Intrusive PWR Algorithm

Step 1: Apply graph decomposition (see Section 3 for details) to identify weakly interacting subsystems in system (1).

Step 2: (Assignment-partitioning process): Partition (1) intomsubsystems (obtained in step 1) leading to system of equations given by

(23). Obtain, λi, λc i and ∑i for each subsystem, i = 1, ..., m. Choose the parameters lsi, lci, Pi.

Step 3: (Initialization of the relaxation process with no coupling effect): Set I = 1 and guess an initial waveform {y0 i (t):t ∈ [0, Ts]} for each

i = 1, ..., m consistent with initial condition (see step 1 in relaxation process described in Section 3.1). Set

(78)

i = 1, ..., m, and solve for yλi,Pi,1(t), t ∈ [0, Ts]} using

(79)

with an initial condition y1i (0) = y0 i (0) on a collocation grid C(o, ni). Determine the gPC expansion yλi,Pi,1i (t, ·) over P-variate

polynomial space Wλi,Pi by computing the expansion coefficients from the quadrature formula (21).

Step 4: (Initialization of the relaxation process, incorporating first level of coupling effect): Set I = 2 and for each i = 1, ..., m, set

(80)

and solve for {y2 i (t), t ∈ [0, Ts]} from

(81)

with an initial condition y2 i (0) = y0i (0), over a collocation grid C(l, ni + nci). Obtain the expansion yi,Pi,2i (t, ·) using (75). From now on

we shall denote the solution vector of the ith subsystem at Ith iteration by yi,Pi,Ii.

Step 5: (Analyzing the decomposed system at the Ith iteration): For each i = 1, ..., m, set

(82)

and solve for {yI (t) : t ∈ [0, Ts]} from

(83)

with initial condition yI i (0) = y0i (0) over a collocation grid C(l, ni + nci). Obtain the expansion yi,Pi,I i (t, ·) using the expansions (75).

Step 6: (Iteration): Set I = I + 1 and go to step 5 until satisfactory convergence has been achieved.

4.5 Convergence of PWR

Below, we give the convergence proof for the iterative PWR approach. The proof is based on showing that the AGP system is Lipschitz if the original systems is Lipschitz (see condition 2), and then invoking the standard WR convergence result (Section 3.1).

Proposition 4.1. Convergence of PWR: The intrusive and nonintrusive PWR algorithms described in Sections 4.3.1 and 4.4.1, respectively, converge.

Proof: We prove the result for intrusive PWR. By construction, because nonintrusive PWR algorithm solves the AGP system (69) in a different way, the convergence result holds for nonintrusive PWR as well. Consider the AGP system (69) and let

(84)

and

(85)

Let for a given k = 1, ..., n, i = 𝒟(k), then

(86)

and 𝒫i(xl,∑,P) = [𝒫i,1(xl,∑,P1), ..., 𝒫i,n(xl,∑,Pn)], for l = 1, 2 and for simplifying notation we have dropped subscripts on P vectors. Then for each k = 1, ..., n, i = 𝒟(k), j = 1, ..., Ni,

(87)

where,

(88)

For a given i = 1, ..., n and j = 1 ..., N𝒟(i) , let Lg ij = [L𝒟(i)𝒟(1)j1 ...L𝒟(i)𝒟(1) jN𝒟(1) ...L𝒟(i)𝒟(n) j1 ...L𝒟(i)𝒟(n) jN𝒟(n)] and define

(89)

then,

(90)

where L = ||Lg|| is a matrix norm of Lg. Hence, the system (69) is Lipschitz. Thus, given the original system (1) is Lipschitz [condition (2)], the approximate system (69) is also Lipschitz as shown above (90). Hence, by proposition 3.1 (adaptation of Theorem 5.3 in [18]), we conclude that PWR converges.

The question of how the decomposition of the system and the choice of the PWR algorithm parameters P, ls, lc influence: (1) the rate of convergence of PWR, and (2) the approximation error [due to the truncation introduced in the AGP system (69), the use of adaptive collocation grid, i.e., condition (77) and computation of the modal coefficients by the quadrature formula], needs to be further investigated.

4.6 Scalability of PWR with Respect to Standard PCM

The scalability of nonintrusive PWR relative to PCM is shown in Fig. 2, where the ratio 𝓡F/𝓡I indicates the computation gain over standard full grid PCM approach applied to the system (1) as a whole. Here 𝓡F = lp is the number of deterministic runs of the complete system (1), which comprises m subsystems each with pi, i = 1, ..., m uncertain parameters, such that p = ∑mi=1 pi and l denotes the level of full grid. Similarly, 𝓡I = 1 + ∑mi=1 lpis + Imax(∑mi=1 lpisji lpic) is the total computational effort with the PWR algorithm (including additional effort incurred due to several iterations), where Imax is the number of PWR iterations. Clearly, the advantage of PWR becomes evident as the number of subsystems m and parameters in the network increases. Note that the linear scaling depicted in Fig. 2 is only for the ratio 𝓡F/𝓡I. Moreover, PWR is inherently parallelizable.

FIG. 2: Scalability of PWR algorithm with respect to standard PCM, when implemented with full grid collocation as the subsystem level UQ method, with pi = 5, ∀i = 1, ..., m, l = ls = 5, lc = 3 and Imax = 10, 50, 100: these numeric values are chosen for illustration purposes. The computational gain of PWR becomes insensitive to Imax (number of iterations), as the number of subsystems m increases.

5. EXAMPLE PROBLEMS

In this section we illustrate intrusive and nonintrusive PWR through several examples of linear and nonlinear networked systems with increasing number of uncertain parameters. While most examples are of ODEs, we also give an example application of PWR to an algebraic system. This illustrates how in principle one can extend the application of PWR to DAEs, just like the WR approach extends to DAEs [15]. Through some examples we study how the strength of interaction between subsystems affects the convergence rate and the approximation error of PWR. In one of the examples related to the building model, we also show how time-varying uncertainty can be incorporated into standard UQ framework by using the Karhunen-Loeve expansion. In all the examples, we compare the solution accuracy of PWR with other UQ approaches (e.g., Monte Carlo and quasi-Monte Carlo methods, QMC), and wherever appropriate mention computation gain offered by PWR over the standard application of gPC and PCM.

5.1 Stability Problem

We first consider a simple system, with two states (x1, x2) ∈ ℝ2,

(91)
(92)

where c, v1, v2 are fixed parameters, and a, b are uncertain with Gaussian distribution 𝒢 and tolerance 20% (i.e. σ = 0.2μ, where μ is the mean and σ is the standard deviation of 𝒢). The parameter c determines the coupling strength between two subsystems described by the two equations. The output of interest is the stability of the system, which is determined by λm, the maximum eigenvalue of the Jacobian

where x10, x20 is the equilibrium point satisfying

(93)

Figure 3 shows the UQ results obtained by applying PWR [with ls = 5, lc = 3, and P1 = (5, 3), P2 = (3, 5)] to iteratively solve the algebraic system (93). We make comparison with the true [to imply a more accurate result obtained by solving the complete system (93)] distribution of λm obtained by using a full collocation grid on the parameter space (a, b) with la = 5, lb = 5, P = (5, 5). PWR converges to the true mean and variance as shown in the left and right panels of Fig. 3 for two different values of c. As the coupling strength c increases (see the right panel in Fig. 3), the number of iterations required for the convergence increases, as expected.

FIG. 3: Left panel: Convergence of mean (μ) and variance ( σ) of λm for c = 0:1. Right panel: Convergence of mean and variance of λm for c = 2:8. Black line indicates the true values.

5.2 Building Example

For energy consumption computation, a building can be represented in terms of a reduced order thermal network model of the form [29]

(94)

where TRn is a vector comprising internal zone air temperatures and internal and external wall temperatures; A[u(t); ξ] is the time-dependent matrix with ξ being parameters, u(t) is control input vector (comprising zone supply flow rate and supply temperature), vectors Qe = [Tamb(t), Qs(t)]T represent the external (outside air temperature and solar radiation), and Qi is the internal (occupant) load disturbances.We consider the problem of computing uncertainty in building energy consumption due to uncertainty in building thermal-related properties and uncertain disturbance loads. These uncertainties can be categorized into: (i) static parametric uncertainty which include parameters such as wall thermal conductivity and thermal capacitance, heat transfer coefficient, window thermal resistance etc.; and (ii) time-varying uncertainties which include the external and internal load disturbances.

Recall that the traditional UQ approaches and PWR, which builds on them, can only deal with parametric uncertainty. To account for time-varying uncertain processes, we employ the Karhunen-Loeve (KL) expansion [30]. The KL expansion allows representation of second-order stochastic processes as a sum of random variables. In this manner, both parametric and time-varying uncertainties can be treated in terms of random variables.We next demonstrate both intrusive and nonintrusive PWR methods.

5.2.1 Two-Zone Example

We first consider a simplified two-zone building model as shown in Fig. 4. Here the state T is a 10-dimensional vector comprising internal wall temperatures and the internal zone air temperatures, where we have assumed that the outer wall surfaces are held at ambient temperature. We also assume that the ambient temperature and solar load are deterministic fixed quantities and there is no internal occupant load. Thus, in computing the uncertainty in the zone temperatures, we only consider parametric uncertainty. Specifically, we assume that the heat transfer coefficient and the thermal conductivity of the walls in each zone have standard deviations of 10% around their nominal values of 3:16 W/m2/K and 4:65 W/m/K, respectively. Thus, locally each zone is affected by two uncertain parameters, with heat transfer coefficient being a common (i.e., same) parameter and thermal conductivity being the other. Using complete Galerkin projection with Pi = (2, 2, 2), i = 1, 2, gives rise to a 60-state ODE model. To apply WR/AWR to this system, we first identify the weakly interacting states. By construction the two zones weakly affect each other, which is identified by the spectral clustering [13] (or wave equation-based clustering [16]) applied to the system (94). This decomposition is imposed on the complete Galerkin system, as explained in Section 4.3.1. As expected, we found that if one applies spectral clustering to the complete Galerkin system instead, one recovers the same decomposition.

FIG. 4: Diagram of the two-zone thermal model of a building. Tamb(t) = 293 K in this example.

Treating 1000 Monte Carlo (MC) samples as the truth, we compare the results of a simulated full Galerkin projected system using both standard waveform relaxation [18] as well as adaptive waveform relaxation [17] in Fig. 5(a). AWR provides a speed-up by a factor of ≈ 12. In Fig. 5(a), one can see that the complete Galerkin projection predicts the same temperature variation over 8 h as MC-based methods.

FIG. 5: (a) Comparison of Monte Carlo, complete Galerkin projection, and approximate Galerkin projection. (b) Normalized error in waveform relaxation as a function of iteration count with increasing coupling. Complete Galerkin and approximate Galerkin are shown. Approximate Galerkin system is found to have greater error as a function of iteration number.

As explained before, one can further exploit the weak interaction between the two zones to reduce the overall number of equations in Galerkin projection. To construct the AGP system, we reduce the order of expansion for the random parameters indirectly affecting each zone so that P1 = (2, 2, 1) and P2 = (1, 2, 2). With this, the number of equations in Galerkin projection reduces from 60 to 50. The resulting solution is shown in Fig. 5(a). We see that the error starts to grow as time increases. However, over 8 h the max error in the room temperatures is 5 × 10−2 K. Thus, despite reducing the computational effort, one can still get a fairly accurate answer.

Figure 5(b) shows the effect of coupling (which is the reciprocal of the coefficient of thermal conductivity of the internal wall) on errors introduced in complete and approximate Galerkin projections. As expected, the approximate Galerkin projection has a higher error [given by ET(t)] than complete Galerkin projection [given by EC(t)]. Moreover this error is more pronounced at low iteration numbers. From the figure, it also clear that as the coupling increases, the number of iterations required for obtaining same solution accuracy increases. For further discussion on the relationship between the coupling and number of iterations, see [17].

5.2.2 Multizone Example

In this section, we consider a larger 6-zone building thermal network model with 68 states. This model admits a decomposition into 23 subsystems, as revealed by the spectral graph approach [see Fig. 6(b)]. This decomposition is consistent with three different timescales (associated with external and internal wall temperature, and internal zone temperatures) present in the system, as shown by the three bands in Fig. 6(a).

FIG. 6: (a) Bands of eigenvalues of time-averaged A(t, ξ) for nominal parameter values. (b) First spectral gap in graph Laplacian revealing 23 subsystems in the network model.

Next we demonstrate the nonintrusive PWR approach to compute uncertainty in energy consumption due to both parametric uncertainty and time varying uncertain loads. As described earlier, we use the KL expansion to transform time-varying uncertainty into parametric form.

KL Expansion [30]: Let {Xt = X(ξ, t), t ∈ [a, b]} be a quadratic mean-square second-order stochastic process with covariance functions R(t, s). If {φn(t)} are the eigenfunctions of the integral operator with kernel R(·, ·) and {λn} the corresponding eigenvalues, i.e.,

(95)

then

(96)

where X(t) is the mean of the process and the limit is taken in the quadratic {an} sense. The random coefficients fang satisfy

(97)

and are uncorrelated E[aman] = δmn. The basis functions also satisfy the orthogonality property

(98)

and the kernel admits an expansion of the form

(99)

Generally, an analytical solution to the eigenvalue problem (95), also known as the Fredholm equation of the second kind, is not available. Several numerical techniques have been proposed; we used the expansion method described in [31].

To apply the KL expansion to the building problem, we assume that the stochastic disturbances [Tamb(t), Qs(t), Qint(t)] are Gaussian processes. This guarantees that the random variables an in the KL expansion are independent Gaussian random variables with a zero mean [31]. Let the joint distribution of a nonstationary Gaussian process be

(100)

where ρ(t, s) is the correlation coefficient and is related to the covariance kernel as

(101)

We assumed the processes Tamb(t), Qs(t) to have a stationary exponential correlation function

(102)

with a constant variance σ2 and a constant correlation timescale Tc. For the internal occupancy load Qint(t) we constructed R(t, s) as follows. For a typical office building, we know that the occupancy load is negligible with low variance during the early and later parts of the day. During peak hours in the middle of the day the occupant load can show significantly high variability. To capture this effect we divided the normalized time domain [0, 1] = [0, t1] ∪ (t1, t2) ∪ (t2, 1) and obtained the desired variation by choosing [in expression (101)]

(103)

with the correlation timescale

(104)

and parameter a controls the slope of the tanh function. Figure 7 shows the covariance kernel for external [Tamb(t), Qs(t)] and internal Qint(t) loads. For the choice of parameters indicated in Fig. 7, we found, using the expansion method [31] with Legendre polynomials as the basis functions, that KL expansion up to order 3 and up to order 6 can capture more than 90% of total variance, for internal and external loads, respectively. Finally, note that these external variables act like global variables which affect all zones in the building model. Such global variables can be easily handled in the PWR framework by treating them as local variables [and hence including them in the local parameter vector defined in (27)] for each subsystem (here the different building zone) they affect.

FIG. 7: (a) Covariance kernel (102 for external load with Tc = 0.1 and σ = 0.1. (b) Covariance kernel (103) for internal load with t1 = t2 = 0.3, a = 20, σ = 0.1.

In UQ computation, we considered the effect of 14 random variables comprising external wall thermal resistance in the six zones, and first dominant random variable obtained in the KL representation of internal load (for each zone) and first two dominant random variable obtained in KL expansion for solar load. Figure 8 show the nonintrusive PWR results on the decomposed network model. As is evident, the iterations converge rapidly in two steps with a distribution close to that obtained from QMC (using a 25,000-sample Sobol sequence) applied to the 68 thermal network model (94) as a whole.

FIG. 8: Histogram of building energy computation for two iterations in PWR. Also shown is the corresponding histogram obtained by QMC for comparison.

5.3 Coupled Oscillators

Finally, we consider a coupled phase only oscillator system which is governed by nonlinear equations

(105)

where n = 80 is the number of oscillators, ωi, i = 1, ..., n is the angular frequency of oscillators, and K = [Kij] is the coupling matrix. The frequencies ωi of every alternative oscillator, i.e., i = 1, 3, ..., 79, is assumed to be uncertain with

a Gaussian distribution with 20% tolerance (i.e., with a total p = 40 uncertain parameters); all the other parameters are assumed to take a fixed mean value. We are interested in the distribution of the synchronization parameters R(t) and phase φ(t) defined by R(t)eφ(t) = 1/NNj=1 eixj(t). Figure 9 shows the topology of the network of oscillators (left panel), along with the eigenvalue spectrum of the graph Laplacian (right panel). The spectral gap at 40 implies 40 weakly interacting subsystems in the network.

FIG. 9: Left panel shows a network of N = 80 phase only oscillators. Right panel shows spectral gap in eigenvalues of normalized graph Laplacian, that reveals that there are 40 weakly interacting subsystems.

Figure 10 shows UQ results obtained by application of nonintrusive PWR to the decomposed system with ls = 5, lc = 2. We make a comparison with QMC, in which the complete system (105) is solved at 25,000 Sobol points [32]. Remarkably, the PWR converges in 4–5 iterations, giving very similar results to that of QMC. It would be infeasible to use full grid collocation for the networks as a whole, since even with lowest level of collocation grid, i.e., l = 2 for each parameter, the number of samples required become RF = 240 = 1:0995e + 012!.

FIG. 10: Convergence of mean of the magnitude R(t) and phase φ(t), and the respective histograms at t = 0:5.

6. CONCLUSION AND FUTURE WORK

In this paper we have proposed an uncertainty quantification approach which exploits the underlying dynamics and structure of the system. Specifically, we considered a class of networked system whose subsystems are dynamically weakly coupled to each other. We showed how these weak interactions can be exploited to overcome the dimensionality curse associated with traditional UQ methods. By integrating graph decomposition and waveform relaxation with generalized polynomial chaos and probabilistic collocation framework, we proposed an iterative UQ approach which we call probabilistic waveform relaxation. We have developed both intrusive and nonintrusive forms of PWR. We proved that this iterative scheme converges under weak assumptions and illustrated it in several examples with promising results. In this work, we have also introduced an approximate Galerkin projection that, based on the results of graph decomposition, computes “strong” and “weak” influence of parameters on states. An appropriate order of expansion, in terms of the parameters, is then selected for the various states.

This approach is shown to accelerate UQ considerably. Several questions need to be further investigated; these include the effect of the choice of parameters on the rate of convergence and approximation error of the PWR algorithm. In order to exploit multiple timescales that may be present in a system, multigrid extensions [33] of PWR may also be desirable. In future work, we also intend to test our approach on much larger dynamic systems that arise in building applications

ACKNOWLEDGMENTS

This work was in part supported by DARPA DSO (Dr. Cindy Daniell PM) under AFOSR contract no. FA9550-07-C- 0024 (Dr. Fariba Fahroo PM). Any opinions, findings and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the AFOSR or DARPA.

REFERENCES

1. Fishman, G., Monte Carlo: Concepts, Algorithms, and Applications, Springer-Verlag, New York, 1996.

2. Niederreter, H., Random Number Generation and Quasi-Monte Carlo Methods, SIAM, Philadelphia, PA, 1992.

3. Sloan, I. H. and Joe, S., Lattice Methods for Multiple Integration, Oxford University Press, Oxford, 1994.

4. Xiu, D. and Karniadakis, G., The Wiener-Askey polynomial chaos for stochastic differential equations, SIAM J. Sci. Comput., 24:619–644, 2002.

5. Xiu, D., Efficient collocational approach for parametric uncertanity analysis, Commun. Comput. Phys., 2(2):293–309, 2007.

6. Hegland, M., Heinrich, S., and Sloan, I. H., Guest Editors′ Preface, J. Complex., 26(5):407, 2010.

7. Gerstner, T. and Griebel, M., Numerical integration using sparse grids, Numer. Algorithms, 18(3-4):209–232, 1998.

8. Foo, J., Wan, X., and Karniadakis, G., The multi-element probablistic collocation method (ME-PCM): Error analysis and applications, J. Comput. Phys., 227(22):9572–9595, 2008.

9. Sobol, I. M., Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates, Math. Comput. Simul., 55:271–280, 2001.

10. Callier, M., Chan, W. S., and Desoer, C. A., Input-output stability theory of interconnected systems using decomposition techniques, IEEE Trans. Circuits Syst., 23(12):714–729, 1976.

11. Georgiou, T. T. and Smith, M. C., Linear systems and robustness -a graph point of view, Lect. Notes Control Inf. Sci., 183:114– 121, 1992.

12. Mezic, I., Coupled nonlinear dynamical systems: Asymptotic behavior and uncertainty propagation, Proc. of 43rd IEEE Conf. on Decision and Control, Atlantis, Paradise Island, Bahamas, Dec. 14–17, 2004.

13. Luxburg, U. V., A tutorial on spectral clustering, Stat. Comput., 17:395–416, 2007.

14. Chung, F., Spectral Graph Theory, Conference Board of the Mathematical Sciences, Washington, 1997.

15. White, J. and Vincentelli, A. S., Relaxation Techniques for the Simulation of VLSI Circuits, Kluwer International Series in Engineering and Computer Science, Boston, 1987.

16. Sahai, T., Speranzon, A., and Banaszuk, A., Hearing the clusters in a graph: A distributed algorithm, Automatica, 48(1):15–24, 2012.

17. Klus, S., Sahai, T., Liu, C., and Dellnitz, M., An efficient algorithm for the parallel solution of high-dimensional differential equations, J. Comput. Appl. Math., 235(9):3053–3062, 2011.

18. Lelarasmee, E., The waveform relaxation method for time domain analysis of large scale integrated circuits: Theory and applications, Tech. Rep. Memorandum No. UCB/ERL M82/40, UC Berkeley, 1982.

19. Thomas, E., A parallel algorithm for simulation of large neural networks, J. Neurosci. Methods, 98:123–134, 2000.

20. Thomas, E., Parallel Algorithms for Biologically Realistic, Large Scale Neural Modeling, Master′s Thesis, Royal Melbourne Institute of Technology, 1998.

21. Rauber, T. and Rünger, G., Parallel implementations of iterated Runge-Kutta methods, Int. J. Supercomput. Appl., 10(1):62– 90, 1996.

22. Roy, S. and Dounavis, A., Longitudinal partitioning based waveform relaxation algorithm for transient analysis of long delay transmission lines, 2011 IEEE MTT-S International Microwave Symposium Digest (MTT), 2011.

23. McKay, W., Monti, A., Santi, E., and Dougal, R., A co-simulation approach for ACSL-based models, Proc. Huntsville Simulation Conf., Huntsville, AL, Oct. 3-4, 2001.

24. Stoer, M. and Wagner, F., A simple min-cut algorithm, J. ACM, 44(4):585–591, 1997.

25. Fiedler, M., Algebraic connectivity of graphs, Czech. Math. J., 23:289–305, 1973.

26. Fiedler, M., A property of eigenvectors of nonnegative symmetric matrices and its application to graph theory, Czech. Math. J., 25(4):619–633, 1975.

27. Golub, G. and Loan, C. V., Matrix Computations, Johns Hopkins University Press, Baltimore, MD, 1996.

28. Kempe, D. and McSherry, F., A decentralized algorithm for spectral analysis, J. Comput. Syst. Sci., 74(1):70–83, 2008.

29. O′Neill, Z., Narayanan, S., and Brahme, R., Model-based thermal load estimation in buildings, Fourth National Conf. of IBPSA-USA, New York, USA, Aug. 11–13, 2010.

30. Stark, H. and Wood, J. W., Probability, Random Processes and Estimation Theory for Engineers, Prentice-Hall, Englewood Cliffs, NJ, 2002.

31. Huang, S. P., Quek, S. T., and Phoon, K. K., Convergence study of the truncated Karhunen-Loeve expansion for simulation of stochastic processes, Int. J. Numer. Meth. Eng., 52:1029–1043, 2001.

32. Joe, S. and Kuo, F. Y., Remark on algorithm 659: Implementing sobol′s quasirandom sequence generator, ACM Trans. Math. Softw., 29:49–57, 2003.

33. Trottenberg, U., Oosterlee, C. W., and Schller, A., Multigrid, Academic Press, London, 2001.

Портал Begell Электронная Бибилиотека e-Книги Журналы Справочники и Сборники статей Коллекции Цены и условия подписки Begell House Контакты Language English 中文 Русский Português German French Spain