Mary R Schurgot
email: mschurgo@stevens.edu

Cristina Comaniciu
email: ccomanic@stevens.edu

Katia Jaffrès-Runser
email: katia.jaffres-runser@irit.fr

Modeling Opportunistic Data Delivery in Dynamic Wireless Networks

  L'archive ouverte pluridisciplinaire

I. INTRODUCTION

Dynamic wireless ad hoc networks are characterized by intermittent connectivity and pairwise contact enabling opportunistic communication. In these ad hoc networks an endto-end path does not exist at any given time, but rather links between nodes change dynamically with varying network conditions, usually based on the mobility in the network. Additional terminology used to describe such networks includes delay tolerant networks, disruption tolerant networks, intermittently connected ad hoc networks, and opportunistic networks [START_REF] Schurgot | Beyond Traditional DTN Routing: Social Networks for Opportunistic Communication[END_REF].

Message delivery in these networks is enabled by the storecarry-forward routing paradigm. Nodes leverage mobility to construct an asynchronous end-to-end path between a source and destination [START_REF] Boldrini | Exploiting users' social relations to forward data in opportunistic networks: The HiBOp solution[END_REF]. When a relay node carrying a message comes within transmission range of another node, the two nodes are said to be in contact. The decision to be made at this time is whether to forward the message to the encountered node or replicate the message by continuing to store the message and transferring a copy to the encountered node. Opportunistic routing protocols employ algorithms which decide to whom to forward and how much to replicate. These routing approaches range from naive flooding to networkaware predictive forwarding.

Flooding protocols operate with the obvious tradeoff that the more you flood the better the chance to reach the destination and in a shorter time, but in turn you waste network resources, in terms of bandwidth, energy and storage. The design of most of the flooding based protocols does not account for the secondary effects of contention for the limited buffer space and the limited network capacity of typical wireless networks, thus resulting in message dropping and retransmissions, with a detrimental effect on performance [START_REF] Lindgren | Probabilistic Routing in Intermittently Connected Networks[END_REF] [START_REF] Tseng | The Broadcast Storm Problem in a Mobile Ad Hoc Network[END_REF].

In contrast, network-aware algorithms aim to predict suitable relays for a destination. Simply adapting to the network conditions to support any opportunistic communication is no longer sufficient; instead, efficient forwarding schemes have been designed to balance the overhead from redundant copies with successful delivery and minimal delay. A new class of opportunistic routing protocols considers the underlying social network of the operators carrying the wireless devices to aid in the identification of efficient carriers for specific destinations. Social-based routing is rooted in the observation that the social connections between users vary more slowly than the wireless links; thus they serve as a good predictor of future contact [START_REF] Hui | BUBBLE Rap: Social-based Forwarding in Delay Tolerant Networks[END_REF]. Another motivation for this type of design is to quantify the wireless network structure to understand the role of nodes in connecting the sparse network for data delivery [START_REF] Schurgot | Beyond Traditional DTN Routing: Social Networks for Opportunistic Communication[END_REF].

As the amount of knowledge available to the protocol increases, network performance with respect to average delay and delivery ratio improves [START_REF] Jain | Routing in a Delay Tolerant Network[END_REF]. With perfect knowledge of connectivity-past, present, and future state informationan optimization problem can be formulated. Solutions for benchmarking opportunistic routing protocols propose various oracle type schemes (e.g. [START_REF] Spyropoulos | Efficient Routing in Intermittently Connected Mobile Networks: The Single-Copy Case[END_REF]), which assume knowledge of all encounters in time and space to construct minimum cost routes. The authors in [START_REF] Spyropoulos | Efficient Routing in Intermittently Connected Mobile Networks: The Single-Copy Case[END_REF] extend their analysis of message forwarding to the multiple copy case with replication [START_REF]Efficient Routing in Intermittently Connected Mobile Networks: The Multiple-Copy Case[END_REF]. However, awareness of the entire wireless ad hoc network at individual nodes is not realistic in practice, let alone knowledge of the future network state.

Our approach to understanding data delivery is to use the statistical properties of node mobility to model data exchange during pairwise contact. A link exists between two nodes with some probability. This probability can be determined from the proposed mobility models in the literature [START_REF] Karamshuk | Human Mobility Models for Opportunistic Networks[END_REF].

Modeling Opportunistic Data Delivery in Dynamic Wireless Networks

Our model draws on the results that the intercontact time between successive node encounters follows a power law distribution with exponential tail [START_REF] Karagiannis | Power Law and Exponential Decay of Intercontact Times between Mobile Devices[END_REF]. Karagiannis et al. note that the exponential tail is an important characteristic based on the time scale of opportunistic networks. This statistical characterization of mobility has also been proposed in [START_REF] Wang | Social Structure Based Routing of Intermittently Connected Network Using Contact Information[END_REF] and [START_REF] Singh | Optimal Forwarding in Delay Tolerant Networks with Multiple Destinations[END_REF] to characterize the network graph. The authors in [START_REF] Wang | Social Structure Based Routing of Intermittently Connected Network Using Contact Information[END_REF] seek to define a social-based routing protocol, which is not a goal here. The authors in [START_REF] Singh | Optimal Forwarding in Delay Tolerant Networks with Multiple Destinations[END_REF] address the problem of distributing a message to multiple destinations where nodes have full knowledge of the number of relays and destinations possessing a copy of the message. In this work, we address a similar problem as in [START_REF] Singh | Optimal Forwarding in Delay Tolerant Networks with Multiple Destinations[END_REF] but we do not make full knowledge assumptions and consider only the case where a single copy of a packet is present in the network.

In this work, our goal is not to propose a new routing protocol for opportunistic communication. Instead, we present a model to capture the exchange of data between nodes forming a dynamic wireless ad hoc network. Quantifying this data exchange is a first step towards a one-to-one protocol performance comparison as the forwarding probabilities can be selected to implement the forwarding policies for any practical routing protocol for opportunistic networks. In future work, the proposed model can be leveraged to formulate an optimization framework to derive better forwarding strategies for routing in intermittently connected networks. The purpose of this paper is to introduce and assess our proposed analytic model for data delivery in opportunistically connected wireless networks. Model results are validated using OMNet++ [START_REF]Omnet: Objective modular network testbed[END_REF] simulations. Model accuracy has been tested for several mobility scenarios for a 4 and 10-node network. Our simulation results show that relative error between the proposed analytical model and simulation is very low (i.e. 0.075 +/-0.04% for the 10-node network), thus validating our model of data flows.

The paper is organized as follows. The network model is presented in Section II. In Section III we derive our proposed analytic model for packet exchange. We present analysis and simulation results in Section IV prior to concluding in Section V.

II. NETWORK MODEL

The network of interest is an intermittently connected mobile network characterized by pairwise contacts between nodes in transmission range. In this model, the network is composed of a set V of N mobile nodes, decomposed into source, destination and relay nodes. Messages are routed in this network from one source node S ∈ V to one destination node D ∈ V in a multi-hop fashion using a store-and-forward paradigm. Source and destination do not participate in data relaying. As nodes come into contact, they exchange data to enable message delivery following the protocol model defined hereafter. In this work, we assume as well that the wireless link between a pair of nodes in contact has an infinite capacity, meaning all packets to be forwarded from one node to the other one are perfectly and instantaneously delivered.

A. Mobility model

D ij i j

Fig. 1. The link between i and j is a function of the intercontact time d.

The opportunistic contacts between two nodes are commonly characterized using their intercontact time and their contact duration statistics. In this model, since we assume infinite transmission capacity when two nodes are in contact, we will not account for the contact duration in our model.

Intercontact time between two nodes i and j is characterized using a random variable D ij (cf. Fig 1) which is exponentially distributed with rate 1/d ij following the results of [START_REF] Karagiannis | Power Law and Exponential Decay of Intercontact Times between Mobile Devices[END_REF]. As such, average intercontact time of i and j is d ij > 0.

B. Protocol model

When a node i encounters another node j, it has two options: do nothing or forward the entire contents of its buffer to j. As such, the forwarding decision of a node i to node j is a binary discrete random variable F ij that follows a Bernouilli distribution of probability x ij . If F ij equals one, the content of the buffer of i is forwarded to j.

In the rest of the paper, we will call x ij the forwarding probability between i and j. This probability is conditional to the fact that i and j meet. In other words, if j is the next node encountered by i and x ij = 0, node i empties its buffer into the buffer of j with probability x ij . Note that if x ji = 0, node j empties its buffer into the buffer of i with probability x ji during the same encounter.

We introduce the following terminology:

• A positive encounter for node i is an encounter where it receives packets from another node. It increases the number of packets stored in its buffer. An encounter is positive for node i if x ji = 0. • A negative encounter for node i is an encounter where it forwards packets to another node. It decreases the number of packets stored in its buffer. An encounter is negative for node i if x ij = 0. Note that an encounter can be both positive and negative.

C. Node model

Nodes can buffer a maximum of K packets. We assume in the rest of this paper that all packets have the same size which is normalized to one.

To characterize the flows of packets transiting through node i, we define the following random variables:

Q in i the
i ).
The destination node is not forwarding any packets. As such, only Q in D is of interest here. We assume that the source node is constantly generating packets with a constant rate R S and is able to empty its buffer at the time of negative encounters. As such, Q in S represents the random variable of the number of packets generated and, similarly to relay nodes, Q out S represents the number of packets forwarded by S at its next negative encounter.

III. ANALYTIC DTN MODEL

The aim of this paper is to propose an analytic model to derive the expected values Qout i for all relays, Qout S and Qin D , as a function of the forwarding probabilities and the mobility of the nodes. As we will highlight in Section IV, forwarding probabilities can model different common DTN routing strategies. In future work, the proposed model will be considered as a building block for deriving optimal routing strategies in an optimization problem where forwarding probabilities are varied. This section introduces first some preliminary calculations and then presents our analytic DTN model.

A. Preliminary calculations a)

Average time until the next encounter of node i: or the minimum time between encounters, is the minimum of a set of exponential random variables, which is also exponentially distributed with a rate equal to the sum of the rates for all 1 dij . It follows that the average time until node i's next encounter is:

∆ i = 1 j =i 1 dij (1)
b) Probability P ij that node i meets node j next: The probability that j is the next node encountered by node i can be calculated based on the properties of the exponential distribution. The probability that node j is the one with the minimum time between encounters is equal to the rate between i and j divided by the sum of all the rates. P ij can thus be defined in terms of ∆ i :

P ij = ∆ i d ij (2)
c) Average time until the next negative encounter of node i: A negative encounter for node i meeting j arises when x ij is positive. As such, the number of forwarding trials from i to j is modeled as a geometric random variable with probability x ij , supported on the set {1, 2, 3, . . . }. The average number of encounters with j node i experiences before a negative encounter happens is 1/x ij . Thus, if these two nodes have an intercontact time d ij , the intercontact time between negative encounters is:

s ij = d ij x ij , ∀ x ij > 0 (3)
The intercontact times between negative encounters retains the properties of the exponential distribution, but is now adjusted with rate 1 sij . The average time until node i's next negative encounter is:

∆ out i = 1 j ={i,S} 1 sij (4)
This value is infinite for nodes that never choose to forward packets. However, in this work, we assume a relay node will always forward to the destination. Thus, ∆ out i exists for all i. d) Average time until the next positive encounter of node i: Similarly, a positive encounter occurs when x ji > 0, i.e. node j's forwarding decision to i is positive and the intercontact time between positive encounters with j is:

u ij = d ij x ji , ∀ x ji > 0 (5)
Thus, the minimum average time until the next positive encounter for node i is:

∆ in i = 1 j ={i,D} 1 uij (6)
This value is infinite for nodes that never receive packets from any other node. For nodes having j ={i,D} x ji > 0, ∆ in i exists. e) Probability P in ij that node i meets node j in a positive encounter: The probability that node j is the next node involved in a positive encounter with i can be calculated as in Eq. [START_REF] Boldrini | Exploiting users' social relations to forward data in opportunistic networks: The HiBOp solution[END_REF].

P in ij = ∆ in i u ij (7)
In the case where ∆ in i does not exist for a node i, P in ij is set to 0 for all j.

B. DTN model derivation

In this paper, we provide a derivation for the infinite buffer case, where each node can store an unlimited number of packets. Future work will also address the finite buffer case.

1) Case of the source node: The source S generates packets with a constant rate R s . When S is isolated (i.e. not in contact with any other node), it accumulates the generated packets in its buffer with rate R s . As such, it stores on average QS packets between negative encounters:

QS = R s ∆ out S (8)
with ∆ out S the average time until the next negative encounter of node S.

When S comes in contact with another node j in a negative encounter, it empties its buffer into j's buffer. At this point in time, Qout S is equal to QS .

Qo u t N -2 P i n i , N -2 Qout i Node i Qi Q o u t S P i n i S Q o u t 1 P in i1
Fig. 2. The average flow of packets coming into node i from N -2 relays and the source S at positive encounters.

2) Relay node relationship:

The average number of packets coming into node i on positive encounters, Qin i , is the weighted sum of the average number of packets the other nodes of the network are forwarding on their negative encounters with i (including the source node). This is illustrated in Fig. 2. In the following, the source has an index S ∈ {1, .., N -1}. For clarity, we extract S from the subsequent summations. For node i, it has probability P in ij to receive packets from j as stated in Eq. [START_REF] Spyropoulos | Efficient Routing in Intermittently Connected Mobile Networks: The Single-Copy Case[END_REF]. Thus:

Qin i = Qout S P in iS + j∈[1,..N ],j ={i,S,D} Qout j P in ij (9)
We can derive Qout i , which is equal at the time of a negative encounter to the average number of packets accumulated in the buffer, as a function of Qin i . This relationship originates in the difference between the inter-arrival times of positive and negative encounters. To derive Qout i , we can simply scale Qin i by the ratio R i of negative to positive encounters.

Qout i = Qin i R i ( 10
)
with R i the ratio of the average time between negative encounters to the average time between positive encounters both defined in Eq. ( 4) and ( 6), respectively:

R i = ∆ out i ∆ in i (11)
Indeed, if more negative encounters than positive ones arise (i.e. ∆ out i < ∆ in i ), more packets are being forwarded to the rest of the network, and Qout i < Qin i . Similarly, if on average node i experiences more often positive encounters than negative ones, it builds up packets in its buffer, and Qout i > Qin i . For the case of infinite time until the next positive encounter, i.e. when ∆ in i does not exist, R i equals 0. Since ∆ out i is never infinite, there is no undefined value of R i . Combining Eq. ( 9) and Eq. ( 10), we get (N -2) equations:

Qout i R i = Qout S P in iS + j∈[1,..N ],j ={i,S,D} Qout j P in ij (12)
Note that for nodes with infinite time before a positive encounter (R i = 0), Q in i =0 and computation of Q out i is not relevant for this case.

3) Case of the destination node: The destination node always experiences positive encounters since it is an absorbing node. As such, ∀i ∈ [1, ..N ], i = {D}, x iD ≡ 1; x Di ≡ 0. It follows that the average amount of packets received Qin D per encounter is then:

Qin D = Qout S P DS + i∈[1,..,N ],i ={S,D} Qout i P Di (13)

4) Global network system:

The unknowns of our model are Qout S , Qin D and Qout i defined for all i ∈ V, i = {D, S}. Qout S is derived using the closed form expression of Eq. ( 8). Once Eq. ( 12) is expanded for each i, the linear system constructed by ( 12) and ( 13) is composed of as many equations as unknowns. In vector form

Q o u t S P i n 2 S D 1 2 S Qo u t S P i n 1 S Q o u
A • Q = b: A =        1 R1 -P in 12 . . . -P in 1,N -2 0 -P in 21 1 R2 . . . -P in 2,N -2 0 . . . . . . . . . . . . . . . -P in N -2,1 -P in N -2,2
. . .

1 RN-2 0 -P D1 . . . . . . -P D,N -2 1        (14a) Q =        Qout 1 Qout 2 . . . Qout N -2 Qin D        (14b) b =        Qout S P in 1S Qout S P in 2S . . . Qout S P in N -2,S Qout S P DS        (14c)
Deriving the values for the average number of packets received at D and forwarded by the relays requires us to solve the aforementioned N -by-N linear system. The necessary condition is for the determinant of A to be non-zero.

a) Illustrative example: For the four node network shown in Fig. 3 with a source S, destination D, and two relay nodes:

Qout 1 R 1 -Qout 2 P in 12 = Qout S P in 1S (15a) -Qout 1 P in 21 + Qout 2 R 2 = Qout S P in 2S (15b) Qin D -Qout 1 P D1 -Qout 2 P D2 = Qout S P DS (15c)
This system can be solved if the determinant of A follows:

1 R 1 * 1 R 2 -P in 21 P in 12 = 0
Since ∀(i, j) ∈ V × V, d ij > 0 and all nodes are forwarding data to the destination (i.e. x iD ≡ 1 for i ∈ V, i = D), the determinant is always positive and can be solved.

IV. NUMERICAL RESULTS

A. Simulations setup

We use MATLAB R2011b to perform calculations and solve for Qout S , Qin D and Qout i , ∀i ∈ V, i = {S, D}. OMNeT++ IDE version 4.1 [START_REF]Omnet: Objective modular network testbed[END_REF] is used as a framework to simulate the opportunistic network and validate our theoretical formulation. The simulation traces the number of packets buffered by each node. As in the previously defined model, packets cannot be discriminated and are transferred without knowing their content. As such, we do not model a forwarding rule only applying to packets that are new, i.e. different from the packets in the encountered node's buffer.

We create n 2 self-messages [START_REF] Varga | OMNeT++ User Manual, Omnet++ version 3[END_REF] to indicate contact between a pair of nodes for all n nodes in the network. The self-messages arrive according to an exponential distribution with rate 1 dij . Communication is possible between a node i and j when a self-message selfMessage ij arrives indicating pairwise contact. Upon the arrival of a self-message, node i and j decide with probability x ij and x ji respectively to transfer the entire content of their buffer to their counterpart. The buffer counters of each node are updated accordingly. There is also a source self-message indicating the arrival of a packet every R S time units. The source's buffer is then increased by one.

The simulation stores for each node i the total number of packets forwarded N f (i) and the total number of negative encounters E n (i). Model values of Qout Each of the first 6 scenarios is run for ten trials with different random seeds. The simulation is allowed to run for 10 6 time units. The forwarding probabilities were modified for some example scenarios related to the four node network in Fig. 3. Uniform and non-uniform mobility scenarios are considered with the source generating packets at a rate R S = 1 packets per time unit. In total, 8 different scenarios are presented for the 4-node network. Their properties are depicted in Table III.

B. Four node network

For the first 6 scenarios, only 2-hop communications are allowed. More specifically, the two relay nodes 1 and 2 cannot exchange packets but only deliver the packets received from S to D (i.e. x 12 = x 21 = 0). Depending on the values of x S1 and x S2 , either 1 or 2 nodes are allowed by the source to carry packets to D. In practice, the forwarding decision of the source x Si may be set by external factors like information that i is a good carrier for D as used by the PRoPHET algorithm [START_REF] Lindgren | Probabilistic Routing in Intermittently Connected Networks[END_REF]; or i is a highly central node as used in part by SimBetTS [START_REF] Daly | Social Network Analysis for Information Flow in Disconnected Delay-Tolerant MANETs[END_REF]; or perhaps i and D belong to the same community like the BUBBLE protocol considers [START_REF] Hui | BUBBLE Rap: Social-based Forwarding in Delay Tolerant Networks[END_REF].

Table III gives the analytical results for the four node network with average intercontact times shown. Here we use the term uniform mobility to refer to the scenarios in which the average intercontact time is the same for all node pairs and non-uniform mobility otherwise. The average percent error between analytical and empirical values for Qout 1 , Qout 2 , Qout S and Qin D is plotted in Fig. 4 for each scenario. Percent error is very low for both non-uniform and uniform mobility scenarios. For uniform mobility, error is below 0.2% and for non-uniform mobility, error is below 0.4%.

In scenarios 7 and 8, communication between nodes 1 and 2 is possible. Scenario 1 is modified in Scenario 7 with x 12 = 0.5 as node 1 empties its buffer into node 2 every two encounters. Thus, a 3-hop path from S to 1 and then to 2 exists. In Scenario 8, the source forwards to node 1 and node 2 at every encounter while x 21 = 0.5. Fig. 5 displays the results for one trial of these three hop example scenarios. Error percent is still low for all nodes of the network.

C. Ten node network

Table IV extends the analysis to a ten node network. In this ten node network, d ij varies between 1 and 10 and is chosen randomly. Forwarding decisions are considered as binary: x ij ∈ {0, 1} and chosen arbitrarily. R S is also changed to 10 packets per time unit. The simulation is run 5 times with 5 different seeds. As shown by these example scenarios, the described model provides a reliable framework for capturing data exchange in opportunistic wireless ad hoc networks. During our efforts to validate the model, we gained insight regarding the expected network performance. For the same mobility in scenarios 1-3 and scenarios 4-6, the same number of packets, 2.715 in the non-uniform case and 3.333 in the uniform case, are received at the destination, despite the variety of forwarding decisions in our example scenarios. With a change in forwarding probability, we show that the buffers of individual nodes accumulate more or less packets between negative encounters depending on the ratio R i . Thus, the performance that is expected to vary is the average delivery delay as set by the mobility of the nodes. In future work, as we limit the capacity of relay buffers, we expect to see a dip in delivery. The current model allows us to identify those nodes that build up a large content store. For example, in scenario 2, the average number of packets leaving node 2 is much greater than that of node 1 when it is selected as the relay in scenario 1.

V. CONCLUSION

In this work we have proposed a novel model for data packet forwarding in intermittently connected networks. A forwarding probability parameter for individual nodes was defined to capture forwarding policies of various practical routing protocols. Our analysis has determined the packet accumulation and data flow mechanisms from node to node with an impact on the computation of the expected final delivery of packets. Our analytical results were validated via simulations for some example networks. Future work will consider also the case of finite buffer capacity for individual nodes with an impact on determining the dropping probabilities at intermediate nodes.

The proposed model provides a foundation for performance comparison of existing opportunistic forwarding protocols.

Fig. 3 .

 3 Fig. 3. Data flow in a 4-node network.

i

  are compared to their empirical counterparts N f (i)/E n (i). For the destination node, the number of received packets N r (D) and the number of positive encounters E p (D) is traced. Qin D is compared to N r (D)/E p (D).

Fig. 4 .

 4 Fig. 4. Percent error results for a 4-node network. Node number 0 is the source and node number 3 is the destination.

10 Forwarding

 10 d 01 = d 02 = d 03 = d 12 = d 13 = d 23 = = 0, x 3i = 0, x i3 = 1 for i = {1, scenario S3, S6: x 01 = 0.5, x 02 = 0

Fig. 5 .

 5 Fig. 5. Percent error results for a 4-node network using 3-hop long paths. Node number 0 is the source and node number 3 is the destination.

  number of packets received by node i at its next

	positive encounter, the number of packets forwarded by node i at its next negative encounter, Definition rate source S generates packets average time between i and j meeting Q out i probability i meets next node j average time until next encounter of node i TABLE I Symbol Rs d ij P ij ∆ i NETWORK CHARACTERISTICS minimum time between encounter source packet rate intercontact time probability of contact At the time of a negative encounter for relay i, all packets in the buffer are forwarded to the encountered node if its forwarding decision is successful. With this specific packet exchange protocol, we always have Q i that equals Q out i at negative encounters. As such, the following model only has to characterize the random variables Q out i for all N -2 relay Q Parameter nodes (and not Q

i the number of packets in the buffer of node i.

TABLE III SIMULATION

 III SCENARIOS AND ANALYTICAL RESULTS FOR AVERAGE INCOMING PACKETS FOR THE FOUR-NODE NETWORK Mobility scenarios Scenarios S1-S3: d 01 =2; d 02 =10; d 03 =100; d 12 =8; d 13 =3; d 23 =40 Scenarios S4-S6:

TABLE IV RESULTS

 IV FOR A 10-NODE NETWORK

		Analytical	Simulation % Error
	Qout S Qout 1 Qout 2 Qout 3 Qout 4 Qout 5 Qout 6 Qout 7 Qout 8 Qin D	3.611 2.045 4.331 3.932 4.611 9.178 4.114 6.673 5.484 4.738	3.609 2.045 4.334 3.929 4.605 9.174 4.120 6.664 5.488 4.736	0.054 0.002 0.064 0.061 0.125 0.050 0.168 0.132 0.065 0.028

The first author is now with LGS Bell Labs Innovations, Florham Park, NJ, USA. Email: maschurg@lgsinnovations.com.


Biswarup Bhattacharyya
email: biswarup.bhattacharyya@ifsttar.fr

Eric Jacquelin
email: eric.jacquelin@univ-lyon1.fr

Denis Brizard
email: denis.brizard@ifsttar.fr

Uncertainty quantification of stochastic impact dynamic oscillator using a proper orthogonal decomposition-polynomial chaos expansion technique

A proper orthogonal decomposition (POD) based polynomial chaos expansion (PCE) is utilized in this paper for the uncertainty quantification (UQ) of an impact dynamic oscillator. The time-dependent non-smooth behavior and the uncertainties are decoupled using the POD approach. The uncertain response domain is reduced using the POD approach and the dominant POD modes are utilized for the UQ of the response quantity. Furthermore, the PCE model is utilized for the propagation of the input uncertainties. Two different cases of impact oscillator are considered namely single impact and multiple impact. The contact between two bodies is modeled by Hertz's law. For both the cases, UQ is performed on the projectile displacement, projectile velocity and contact force. A highly non-smooth behavior is noticed for the contact force. For that reason, most number of POD modes are required to assess UQ of contact force. All the results are compared with the Monte Carlo simulation (MCS) and time domain PCE results. Very good accuracies are observed for the PCE and the POD-PCE predicted results using much less number of model evaluations as compared to MCS. As the PCE coefficients are dependent on time, the PCE model is computed at each time-step. On the contrary, for the POD-PCE model, the PCE coefficients are computed for the number of POD modes only: it is much less than the PCE model.

Introduction

Low velocity impact dynamic behavior has been studied for beam structures [START_REF] Pashah | Prediction of structural response for low velocity impact[END_REF][START_REF] Apetre | Low-velocity impact response of sandwich beams with functionally graded core[END_REF] and composite structures [START_REF] Abrate | Modeling of impacts on composite structures[END_REF][START_REF] Wu | Low velocity impact damage in composite sandwich beams[END_REF]. Numerical simulation of these structures using finite element model often leads to a high computational cost and a large time consumption. For that reason, the impact dynamical systems are often modeled as a spring-mass impact oscillator [START_REF] Pashah | Structural response of impacted structure described through anti-oscillators[END_REF][START_REF] Pashah | Prediction of structural response for low velocity impact[END_REF][START_REF] Abrate | Modeling of impacts on composite structures[END_REF] and the contact behavior is modeled using the Hertz's law [START_REF] Stronge | Impact mechanics[END_REF][START_REF] Abrate | Modeling of impacts on composite structures[END_REF]. The variation in the contact force and in the deflection can be noticed for a structure as the projectile speed varies. Due to the workmanship and the fabrication procedure, the variation can also be noticed in the material properties. In a nutshell, the input design parameters should be considered uncertain while designing an impact dynamic oscillator. Hence, uncertainty quantification (UQ) of the response quantities due to such uncertain input parameters is the main objective of this study.

A large number of methods have been developed in the literature for the UQ of stochastic dynamical systems [START_REF] Lucor | Adaptive Generalized Polynomial Chaos for Nonlinear Random Oscillators[END_REF][START_REF] Maitre | Asynchronous Time Integration for Polynomial Chaos Expansion of Uncertain Periodic Dynamics[END_REF][START_REF] Gerritsma | Time-dependent generalized polynomial chaos[END_REF][START_REF] Spiridonakos | Metamodeling of dynamic nonlinear structural systems through polynomial chaos NARX models[END_REF]. The main difficulty with the impact dynamic oscillator is the non-smooth dynamic behavior of the response quantities [START_REF] Pashah | Prediction of structural response for low velocity impact[END_REF]. Along with this, the uncertain behavior of the responses is very difficult to capture using the available methods. For the stochastic dynamical systems, polynomial chaos expansion (PCE) has gained popularity for assessing the stochastic dynamical responses [START_REF] Wan | An adaptive multi-element generalized polynomial chaos method for stochastic differential equations[END_REF][START_REF] Sarrouy | Stochastic analysis of the eigenvalue problem for mechanical 16 Bhattacharyya, VIB-19-1536 systems using polynomial chaos expansion-application to a finite element rotor[END_REF][START_REF] Nechak | Prediction of random self friction-induced vibrations in uncertain dry friction systems using a multi-element generalized polynomial chaos approach[END_REF][START_REF] Bhusal | Uncertainty Quantification Using Generalized Polynomial Chaos Expansion for Nonlinear Dynamical Systems With Mixed State and Parameter Uncertainties[END_REF]. One of the main drawbacks of PCE used in time domain is that the PCE coefficients must be constructed at each time-step for a dynamical systems. For that reason, PCE has been modified in several papers [START_REF] Lucor | Adaptive Generalized Polynomial Chaos for Nonlinear Random Oscillators[END_REF][START_REF] Maitre | Asynchronous Time Integration for Polynomial Chaos Expansion of Uncertain Periodic Dynamics[END_REF][START_REF] Gerritsma | Time-dependent generalized polynomial chaos[END_REF] to properly capture the time-dependent stochastic behavior. In [START_REF] Gerritsma | Time-dependent generalized polynomial chaos[END_REF], the polynomial basis is updated using the stochastic response from the previous time-step because the initial basis function is not optimal with the time. A similar procedure has been adopted in [START_REF] Luchtenburg | Long-time uncertainty propagation using generalized polynomial chaos and flow map composition[END_REF]. To capture the nonlinear dynamical behavior, a different approach has been adopted in [START_REF] Mai | Surrogate models for oscillatory systems using sparse polynomial chaos expansions and stochastic time warping[END_REF] by rescaling the whole time domain to a small time domain and uncertainties were propagated using a principal component analysis (PCA) based PCE. The nonlinear autoregressive with exogenous input (NARX) model along with PCE has been used in [START_REF] Spiridonakos | Metamodeling of dynamic nonlinear structural systems through polynomial chaos NARX models[END_REF] for the UQ of nonlinear dynamical systems. The NARX model has also been used along with the Kriging model in [START_REF] Bhattacharyya | A Kriging-NARX model for uncertainty quantification of nonlinear stochastic dynamical systems in time domain[END_REF].

Although time-dependent nonlinear behavior and randomness were decoupled in the sparse NARX-PCE model [START_REF] Spiridonakos | Metamodeling of dynamic nonlinear structural systems through polynomial chaos NARX models[END_REF], the NARX model requires the knowledge of the physical system. Model order reduction (MOR) techniques also decouple the time-dependent behavior and the randomness. MOR technique includes proper generalized decomposition (PGD) [START_REF] Chinesta | A Short Review on Model Order Reduction Based on Proper Generalized Decomposition[END_REF], proper orthogonal decomposition (POD, also known as PCA) [START_REF] Chatterjee | An introduction to the proper orthogonal decomposition[END_REF][START_REF] Kunisch | Galerkin Proper Orthogonal Decomposition Methods for a General Equation in Fluid Dynamics[END_REF] and dynamic mode decomposition [START_REF] Schmid | Dynamic mode decomposition of numerical and experimental data[END_REF]. PCA based PCE has been used in [START_REF] Blatman | Sparse polynomial chaos expansions of vector-valued response quantities[END_REF] for reducing the dimensionality of the time-dependent stochastic diffusion problem. POD has been used in dimensionality reduction of ordinary differential equation [START_REF] Moore | Principal component analysis in nonlinear systems: Preliminary results[END_REF] and partial differential equation [START_REF] Kunisch | Control of the Burgers Equation by a Reduced-Order Approach Using Proper Orthogonal Decomposition[END_REF][START_REF] Kunisch | Galerkin Proper Orthogonal Decomposition Methods for a General Equation in Fluid Dynamics[END_REF]. POD has been used in [START_REF] Higdon | Computer Model Calibration Using High-Dimensional Output[END_REF] for the dimensionality reduction of high-dimensional outputs. POD has the capability to decompose the time-dependent stochastic behavior using much less number of reduced bases. For that reason, POD is used in this paper in accordance with the PCE model [START_REF] Xiu | The Wiener-Askey polynomial chaos for stochastic differential equation[END_REF] for the UQ of impact dynamic oscillators. This method has already been investigated for the smooth random dynamical systems in [START_REF] Jacquelin | Random dynamical system in time domain: A POD-PC model[END_REF]. In this paper, the non-smooth dynamical behavior is investigated using the POD-PCE surrogate model.

This paper is organized in the following way. The PCE model and the POD approach are discussed in section 2 and section 3 respectively. The combined approach is discussed in section 4. UQ of impact dynamic oscillator is performed in section 5 and the conclusions of this study are discussed in section 6.

Polynomial chaos expansion

A brief overview of the generalized PCE model is discussed here in the context of the dynamical systems. Consider a dynamical system having d independent random variables Ξ ∈ R d ( Ξ j represents the j-th variable). Hence, the time (t) dependent response Y Ξ,t can be represented by the PCE model [START_REF] Xiu | The Wiener-Askey polynomial chaos for stochastic differential equation[END_REF]:

Y Ξ,t = ∞ ∑ i=1 y i (t) Ψ (i) Ξ (1)
where Ψ (i) (•) are the multivariate orthogonal polynomial basis functions and y i (t) are the time-dependent coefficients. In the following, 't' may indicate either the continuous time or the discretized time vector with n t elements. The multivariate orthogonal polynomial basis functions are constructed by the product of the univariate polynomial basis functions as follows:

Ψ (i) Ξ = d ∏ j=1 ψ (i, j) Ξ j (2)
where ψ (i, j) Ξ j is the i-th univariate orthogonal polynomial basis function for the j-th input random variable Ξ j . The univariate orthogonal polynomial bases are directly related to the type of the distributions of the input random variables Ξ j [START_REF] Xiu | The Wiener-Askey polynomial chaos for stochastic differential equation[END_REF] e.g. Legendre polynomial is used for an uniformly distributed random variable and Hermite polynomial is used for a normally distributed random variable. The multivariate orthogonal polynomial bases satisfy the following condition of orthogonality:

Ψ (i) Ψ ( j) = R Ξ Ψ (i) Ξ Ψ ( j) Ξ f Ξ Ξ d Ξ = h 2 δ i j (3) 2 Bhattacharyya, VIB-19-1536
where δ is the Kronecker delta and h is a constant which is 1 for the orthonormal polynomial basis function. f Ξ Ξ is the joint probability distribution function (PDF) of the input random variables.

For the practical implementation of the PCE model, the expression in Equation 1 is truncated with some finite number of terms in the polynomial bases. Hence, the truncated PCE is represented by:

Y Ξ,t ≈ n ∑ i=1 y i (t) Ψ (i) Ξ (4)
where n is the finite number of polynomial bases in the PCE model which can be found as follows:

n = d + P P = (d + P)! d!P! (5)
where 'P' is the maximum degree of the multivariate orthogonal polynomial bases. PCE is utilized in this paper in a non-intrusive way. For the non-intrusive approach, initially the responses are computed as though they were deterministic at some predefined sample points (experimental design points). Therefore, the samples of the input random variables Ξ is given by Ξ ∈ R N×d matrix (N is the number of sample points) and the truncated PCE model can be identified from:

Y (Ξ,t) ≈ n ∑ i=1 Ψ (i) (Ξ) y i (t) (6)
In Equation 6, Y = {Y 1 ,Y 2 , . . . ,Y N } T ∈ R N×n t is a matrix having the responses at the initial N samples and n t is the total number of discretized points in the time domain. Ξ is a matrix having N samples of the d-dimensional input random variables (Ξ i ∈ R 1×d is the i-th sample of the d random variables). Ψ (Ξ) ∈ R N×n is the multivariate orthogonal polynomial basis matrix: element (i, j) is Ψ i, j (Ξ) = Ψ j (Ξ i ). y (t) = {y 1 (t) , y 2 (t) , . . . , y n (t)} T ∈ R n×n t is the time-dependent coefficient matrix of the PCE. One of the most important tasks for the PCE model is the computation of the PCE coefficients. The usual way for computing the PCE coefficients is the ordinary least square (OLS) approach [START_REF] Bhattacharyya | A Critical Appraisal of Design of Experiments for Uncertainty Quantification[END_REF]. The PCE coefficients can be computed using the OLS approach as:

∀t k , k ∈ {1, . . . , n t }; y (t k ) = Ψ (Ξ) T Ψ (Ξ) -1 Ψ (Ξ) T Y (Ξ,t k ) (7)
Hence, for an impact dynamic oscillator, the coefficients must be computed at n t number of time-steps which is computationally expensive. For that reason, POD is used in this paper to decouple the time domain and the randomness.

Proper orthogonal decomposition

Proper orthogonal decomposition is mainly used for the identification of low order bases for a given set of snapshots [START_REF] Sirovich | Turbulence and the dynamics of coherent structures Part I: Coherent structures[END_REF]. The orthogonal bases are selected corresponding to the dominating singular values of the snapshots. In our case, we have a similar problem of time-dependent stochastic response quantities. Recall the responses of a dynamical system for N initial samples Y (Ξ,t) ∈ R N×n t : each column of matrix Y contains the stochastic response at each time-step. Hence, a kind of correlation matrix for the stochastic response quantity [START_REF] Berkooz | The Proper Orthogonal Decomposition in the Analysis of Turbulent Flows[END_REF] can be defined as:

C = Y T Y (8)
where the superscript 'T ' denotes the transpose of the preceding matrix. The orthogonal bases are found from the eigenvalue decomposition of the correlation matrix C which should satisfy the following condition:

CΦ i = λ i Φ i ; i = 1, . . . , n t (9)
where Φ i are the orthogonal bases of the POD, which are also called POD modes. λ i are the eigenvalues of the decomposition and represent a kind of energy [START_REF] Chatterjee | An introduction to the proper orthogonal decomposition[END_REF] related to the corresponding POD orthogonal modes. The POD modes are sorted 3 Bhattacharyya, VIB-19-1536

according to the decreasing magnitude of the corresponding eigenvalues λ i . The POD modes are normalized here as follows:

φ i = Φ i Φ i ( 10
)
where • is the 2-norm of the corresponding vector. The orthogonal POD modes can also be computed using the singular value decomposition (SVD) [START_REF] Chatterjee | An introduction to the proper orthogonal decomposition[END_REF].

Having the orthonormal POD modes, the uncertain response series can be represented by:

Y Ξ,t = n t ∑ i=1 b i Ξ φ T i (t) (11)
where b i Ξ are the random coefficient vectors corresponding to the POD orthogonal modes, which need to be evaluated.

It should be noted from Equation 11that the responses are projected on n t number of orthogonal POD modes. However, in our case N < n t . Hence, the total energy is recovered only in the first N eigenvalues and the responses should be projected using only the first N POD modes. Furthermore, it has been seen in the literature [START_REF] Sirovich | Turbulence and the dynamics of coherent structures Part I: Coherent structures[END_REF][START_REF] Azeez | Proper orthogonal decomposition (POD) of a class of vibro-impact oscillations[END_REF] that only n b < N POD modes are usually required to keep 99.99% of the POD energy. Therefore, using the reduced number of orthonormal POD modes, the random responses can be represented by:

Y Ξ,t = n b ∑ i=1 b i Ξ φ T i (t) (12)

POD-PCE model

The coefficients in Equation 12 are random and are estimated using the PCE model: the random coefficients can be expressed using the truncated PCE model (refer Equation 6) as:

b i Ξ ≈ n ∑ j=1 a i j Ψ ( j) Ξ ( 13
)
Coefficients a i j are identified from the known responses of the N sample points

Y (Ξ,t) = n b ∑ i=1 b i (Ξ) φ T i (t) (14)
where b i (Ξ) ∈ R N×1 is the i-th POD coefficient vector. Each coefficient b i (Ξ) is known and is computed by projecting the random response matrix Y on the corresponding orthonormal POD mode:

b i (Ξ) = Y (Ξ,t) φ i (t) ; i = 1, . . . , n b (15)
Having the random coefficients, the coefficient vector a i j can be computed using the OLS approach (refer Equation 7). Once the truncated PCE coefficients are computed, the combined POD-PCE model is given by:

Y Ξ,t ≈ n b ∑ i=1 n ∑ j=1 a i j Ψ ( j) Ξ φ i (t) (16)

UQ of impact dynamic oscillator

In this section, the UQ of a 2 degree of freedom (DOF) impact dynamic oscillator is performed. The reference solution was computed by the MCS approach using Latin hypercube sample (LHS) points. The full PCE (referred as PCE in the [START_REF] Wu | Low velocity impact damage in composite sandwich beams[END_REF] Bhattacharyya, VIB-19-1536

following) and POD-PCE models were identified to predict the stochastic responses using much less samples than the MCS. The accuracy of the predicted result is illustrated using the mean relative error:

ε = 1 N MCS N MCS ∑ k=1 ε k (17)
where N MCS is the number of samples for the MCS. ε k is the relative error of the k-th sample point which was computed as [START_REF] Mai | Surrogate modeling for stochastic dynamical systems by combining nonlinear autoregressive with exogenous input models and polynomial chaos expansions[END_REF]:

ε k = ∑ n t i=1 Y MCS (t i , Ξ k ) -Y pred (t i , Ξ k ) 2 ∑ n t i=1 [Y MCS (t i , Ξ k ) -ȲMCS (Ξ k )] 2 (18)
where ȲMCS (Ξ k ) is the mean of the MCS computed response series at the k-th sample point. Y pred is the predicted response quantity using a surrogate model.

Problem definition

The 2-DOF nonlinear impact oscillator considered here is shown in Figure 1. The structure is denoted by the mass m st and the projectile is denoted by the mass m p . The projectile is not attached to the nonlinear spring having stiffness k c . The governing differential equation of the 2-DOF impact oscillator is given by:

m st ẍst + c st ẋst + k st x st + f c = 0 m p ẍp -f c = 0 ( 19
)
where the structure mass m st is attached to the ground by a linear spring having stiffness k st and by a damper with viscous damping c st . x st and x p are the displacements of the structure and the projectile, respectively. f c is the contact force between the projectile and the nonlinear spring, and was modeled using the Hertz law [START_REF] Johnson | Contact Mechanics[END_REF][START_REF] Buezas | Propagation of uncertainties and multimodality in the impact problem of two elastic bodies[END_REF]. f c is defined as:

f c = k c (x st -x p ) 3 2 ; x st ≥ x p = 0 ; x st < x p (20)
The initial conditions for the impact oscillator are given by:

x st (0) = 0 ẋst (0) = 0 x p (0) = 0 ẋp (0) = v 0 ( 21
)
where v 0 is the initial velocity of the projectile mass m p . In this paper, m p , k c and v 0 are considered uniformly distributed independent random variables with coefficient of variation 0.1. The parameters of the distribution for all the random variables are given in Table 1. The structure mass m st = 60g and the damping ratio ζ st = 0.5% are considered for this oscillator.

UQ is performed considering the projectile displacement, the projectile velocity and the contact force as the response quantities separately. The behavior of the responses is highly dependent on the number of contacts between the structure and the projectile. For that reason, the number of contacts is tuned here varying k st . In this paper, two cases are considered: (i) single impact (k st = 240MN m -1 ), (ii) multiple impact (k st = 2.4MN m -1 ). The non-smooth dynamical behavior is expected to be higher in case of the multiple impact oscillator.

Case 1: Single impact oscillator

The impact oscillator as introduced in the previous section is tuned here such that the contact between the structure and the projectile occurs only once in the total time period. A high value of k st = 240MN m -1 was considered. The time domain [START_REF] Pashah | Structural response of impacted structure described through anti-oscillators[END_REF] Bhattacharyya, VIB-19-1536

+ 1600 √ 3 MN m -3/2 v 0 -10 + √ 3 -10 - √ 3 m s -1
for this system was decided as t ∈ [0ms, 1ms]. The time integration was performed using the MATLAB function ode45 at a time-step of ∆t = 1µs. As a result, n t = 1001 can be found in the full time domain.

For the reference solution, the stochastic responses were computed by MCS using N MCS = 10 4 model evaluations. On the other hand, the POD-PCE and the PCE models were constructed using much less number of model evaluations. Legendre polynomials were used for the PCE models with degree P = 2 and P = 3. As a result, n = 10 and n = 20 terms were found in the truncated PCE model for P = 2 and P = 3, respectively.

First of all, the influence of N on the accuracy was studied for both the surrogate models. The evaluation of the mean relative errors for all the response quantities is shown in Figure 2. It is seen clearly that with the increase of N the mean relative error reduces for all the response quantities. It is noticeable that for a particular polynomial degree, ε by the PCE and the POD-PCE models are almost the same with N. However, slightly high accuracy is noticed by the POD-PCE model for the contact force with P = 3 and the accuracy of PCE models for the projectile velocity is greater than the POD-PCE models. An accurate identification of a model with n coefficients is generally possible only if N ≥ n. Accordingly, the responses with P = 2 are in good agreement with the MCS responses with N ≥ 10 because n = 10 terms were found in the PCE polynomials, whereas P = 3 requires at least N = 20 model evaluation to compute the PCE coefficients accurately (n = 20). With N = 20 and P = 3, the surrogate model has almost the same accuracy as that with P = 2 and N = 20. Moreover, the increment of the accuracy beyond N = 20 samples is much less for all the response quantities. Therefore, a good accuracy has achieved using N = 20 for all responses. For that reason, all the stochastic responses are presented in the following text using N = 20 only for the single impact oscillator.

The time-dependent mean (µ) and standard deviation (σ) of the contact force are plotted in Figure 3a and Figure 3b, respectively. It is seen that the contact between the structure and the projectile occurs only once during the total time period. The mean and the standard deviation are computed quite well by both the surrogates using N = 20. However, it is to be 6

Bhattacharyya, VIB-19-1536 (c) ε for projectile velocity Fig. 2: Evaluation of the mean relative errors of all the stochastic response quantities for the single impact oscillator noted that the PCE coefficient vector was computed at each time-step for the PCE model (refer Equation 6). Conversely, the time-dependent stochastic contact force was obtained with only n b = 5 POD modes by the POD-PCE model. Hence, the PCE coefficient vector was computed only 5 times for the POD-PCE model which is much less than 1001. Indeed the total number of computed coefficients (n tot = n × n t for PCE and n tot = n × n b for POD-PCE) is shown in Table 2 for both the surrogates, which is 200 times more for the PCE model as compared to the POD-PCE model. The time-dependent errors for the statistical moments were computed using the following expression:

ε (t) = Mom (Y MCS ) -Mom Y pred max |Mom (Y MCS )| (22)
where Mom (•) is the any statistical moment (µ or σ) of the corresponding response quantity. The time-dependent errors for µ f c and σ f c of the contact force are shown in Figure 4a and Figure 4b, respectively: the errors are presented only during the contact and the error is zero beyond the contact. It is seen clearly that the magnitude of error in the prediction of σ f c is higher than that in the prediction of µ f c . Moreover, the POD-PCE model has better predicted σ as compared to the PCE model using P = 3. The overall accuracies of all the models were computed using Equation 17 and are listed in Table 2.

It is noticeable that the errors predicted by the POD-PCE models are slightly less than the PCE models. Moreover, a good accuracy has been achieved by both the surrogate models using much less number of model evaluations, but the POD-PCE model requires only 5 POD modes to achieve a good accuracy. The time-dependent statistical moments of x p are shown in Figure 3c and Figure 3d using N = 20. Good agreement of surrogate predicted results with the MCS results is noticed for both the cases. As the PCE coefficients are time-dependent, the PCE coefficient vector was computed 1001 times. However, 99.99% of the energy was recovered using only 2 POD modes for [START_REF] Lucor | Adaptive Generalized Polynomial Chaos for Nonlinear Random Oscillators[END_REF] Bhattacharyya, VIB-19-1536 Fig. 4: Time-dependent errors for the statistical moments of the stochastic responses of single impact oscillator surrogates is measured using the mean relative error which is listed in Table 2. It is noticeable from the table that almost the same accuracies have been achieved using both the surrogates: for the POD-PCE model, it requires only 2 POD modes to achieve such a high level of accuracy. In a similar way to the previous two responses, the stochastic projectile velocity was also computed using the surrogate models. All the parameters remain the same for the projectile velocity. Time-dependent mean and standard deviation of ẋp [START_REF] Gerritsma | Time-dependent generalized polynomial chaos[END_REF] Bhattacharyya, VIB-19-1536 3e and Figure 3f, respectively. It is seen that the mean velocity is gradually increasing from negative magnitude to positive magnitude during the contact period and it remains constant after the contact. The prediction accuracies of all the surrogates are very high and are very close to each other. For the POD-PCE model, the 99.99% energy level is achieved using only n b = 3 POD modes: indeed only 3 PCE coefficient vectors were computed for the POD-PCE model instead of 1001 PCE coefficient vectors for the PCE model. The time-dependent errors were also computed for the predicted stochastic velocity ẋp and are shown in Figure 4e and Figure 4f: the errors are plotted for the contact period only and beyond this region, the errors are constant. It is noticed that the magnitude of error is higher for σ ẋp as compared to µ ẋp . However, the error for σ ẋp by the POD-PCE model with P = 3 is much less as compared to the others. The overall accuracies of all the surrogates are listed in Table 2. The mean relative errors by all the surrogates are close to each other.

The results obtained for 5 samples are plotted in Figure 5, as well as the results calculated from the identified surrogate models. It is seen that all the responses are predicted quite well by both the surrogate models for all the 5 samples.

The efficiency of all the approaches was measured with the number of sample points and CPU time which are listed in Table 3. The simulation was performed on Intel Core i7 CPU running at 2.90 GHz with 32 GB installed memory. The CPU times presented in Table 3 were averaged for 10 simulations. The surrogate models with P = 3 are presented only. It is seen that the number of sample points for the surrogate models is quite small as compared to the MCS approach. Further, the CPU time for the PCE model is almost 830 times smaller than the MCS approach and almost 21000 times smaller for the POD-PCE model. For the POD-PCE models, the POD was performed by the SVD approach and it depends on the minimum between the number of model evaluations (N) and the number of time-steps (n t ). For the present case, N is quite small as compared to n t . Therefore, the POD approach took approximately 0.001 s to obtain the POD modes. As a result, the efficiency of the POD-PCE model was almost 30 times higher than the PCE model.

Case 2: Multiple impact oscillator

The same oscillator as defined in section 5.1 is considered here with k st = 2.4MN m -1 such that multiple contacts occur between the structure and the projectile. The time domain considered here is t ∈ [0ms, 3ms] and the time integration was performed at a time-step of ∆t = 1µs. Consequently, n t = 3001 number of time-steps can be found in the total time period.

Stochastic responses were computed considering the three input random variables as mentioned in Table 1. For the multiple impact oscillator, MCS was performed using 10 4 model evaluations on LHS points for all the stochastic response quantities. In a similar way to the single impact oscillator, the influence of the number of model evaluations was investigated for the surrogate models. As this system is highly non-smooth and nonlinear in nature, initial sample for the surrogate model was decided as N = 30. The relative error with respect to N for all the response quantities is shown in Figure 6. It is noticeable that the accuracy is increased with the increase of N for all the response quantities. The magnitude of error [START_REF] Spiridonakos | Metamodeling of dynamic nonlinear structural systems through polynomial chaos NARX models[END_REF] Bhattacharyya, VIB-19-1536 for the contact force is quite higher as compared to the other response quantities due to highly non-smooth and nonlinear in behavior. Moreover, a quite good accuracy has achieved using N = 50 samples for all the response quantities. Hence, all the results are presented in following discussion using N = 50 samples. The stochastic behaviors of the contact force are plotted in Figure 7a and Figure 7b: the results are presented magnifying on the x axis. The surrogate models were constructed using N = 50 LHS points and predicted at N MCS = 10 4 samples. It is seen clearly from the mean force that five contacts occur during the whole time period and the maximum contact force occurs during the second contact. Here we can see the much high non-smooth dynamic behavior of the oscillator as compared to the previous one. Despite the much high non-smooth behavior, the mean force was predicted quite well by both the surrogates. However, tiny discrepancies are noticed in the standard deviation of f c . The time-dependent errors were also computed for both the surrogates, which are shown in Figure 8a and Figure 8b for the contact region only. The magnitude of error for σ f c is high as compared to µ f c . Both the surrogates have predicted lower error with P = 3 as compared to P = 2. The level of overall accuracies are listed in Table 4. It is noticeable that a large number of POD modes (n b = 31) was required to keep 99.99% of the POD energy with the retained POD modes due to much high non-smooth nonlinear behavior: the PCE coefficient vector was therefore computed 31 times. However, it is much less than the computation of the coefficient vector for the PCE model: it is 3001 for the PCE model. Indeed n tot for the POD-PCE model is approximately 1% of the time-dependent PCE model.

The stochastic projectile displacement x p was computed by both the surrogates using 50 LHS points. Time-dependent moments of x p are shown in Figure 7c and Figure 7d. Both the moments are predicted quite well using both the surrogates. The time-dependent errors were computed for the moments of x p and are plotted in Figure 8c and Figure 8d. The magnitude of error for σ x p is higher as compared to µ x p . Slightly lower errors were predicted by both the surrogates with P = 3. Along with this, the overall accuracies of all the surrogates are listed in Table 4 exactly the same mean relative errors, however, the PCE coefficient vector was computed 3 times instead of 3001 for the POD-PCE model: it is the main advantage of the POD-PCE model. UQ was performed on the projectile velocity also. The stochastic response behavior in time domain is shown in Figure 7e and Figure 7f: the standard deviation is shown only for the contact duration and the standard deviation is constant beyond this region. A highly non-smooth behavior is noticed in the stochastic behavior of ẋp . The peaks in the standard deviation occurred during the contacts. Such highly non-smooth moments were predicted quite well by both the surrogates using much less model evaluations as compared to MCS. The time-dependent errors for both the statistical moments are plotted in Figure 8e and Figure 8f: the errors are plotted in the contact region and the error is constant beyond this region. It is seen [START_REF] Sarrouy | Stochastic analysis of the eigenvalue problem for mechanical 16 Bhattacharyya, VIB-19-1536 systems using polynomial chaos expansion-application to a finite element rotor[END_REF] Bhattacharyya, VIB-19-1536 Fig. 7: Time-dependent statistical moments of the stochastic responses for the multiple impact oscillator: MCS was performed using N MCS = 10 4 model evaluations that the magnitude of error is significantly higher for σ ẋp as compared to µ ẋp . Low errors are noticed using P = 3 for both the surrogate models during the contacts. The overall accuracies are given in In a similar way to the single impact oscillator, the predicted stochastic responses by the surrogates were plotted and compared with the MCS results at the same 5 samples. The comparison is shown in Figure 9. It seems that the contact force is predicted quite well during the first two contacts, while it is very difficult to predict the contact forces towards the end of the contacts. ε for these 5 samples was found to be the same as the ε for the 10 4 samples. It is noticeable that some negative forces are predicted towards end contacts which are not physically correct. To mitigate this issue, we tried [START_REF] Bhusal | Uncertainty Quantification Using Generalized Polynomial Chaos Expansion for Nonlinear Dynamical Systems With Mixed State and Parameter Uncertainties[END_REF] Bhattacharyya, VIB-19-1536

to identify the POD-PCE model and the PCE model by increasing the samples and the polynomial degree to N = 1000 and P = 15, respectively. Although it was possible to reduce the magnitude of the negative forces to some extent using such large number of samples, some negative forces were always predicted using the surrogate models. The negative forces are not a consequence of the POD but of the PCE as the PCE predicted results follow the POD-PCE predicted results. Indeed predicting a non-smooth function (i.e. contact force) by a smooth function (i.e. PCE model) is always a challenging task. On the other hand, the projectile displacement and the projectile velocity are predicted quite well at the 5 samples using both the approaches due to less non-smooth behavior than the contact force. The accuracies of all the surrogate models have already been shown in Table 4. The efficiency metrics for the responses of the multiple impact oscillator is presented in Table 5. The CPU times were measured by averaging on 10 simulations. It is noticeable that the required CPU times for both surrogate models are quite low as compared to the MCS. The number of sample points is quite small as compared to the number of time-steps. As a result, the required CPU time to obtain the POD modes was approximately 0.008 s. The CPU times for the POD-PCE models are almost 16 times and 10000 times lower than that for the PCE models and for the MCS approach, respectively. Therefore, the POD-PCE model is quite efficient as compared to the other approaches. A stochastic impact oscillator has been investigated in this paper by POD-PCE model. The contact between the structure and the projectile was modeled using the Hertz's law. The responses of projectile were considered for UQ. The timedependent UQ can be performed using the conventional PCE surrogate model. However, the PCE coefficient vector has to be computed at each time-step. For that reason, the time-dependent non-smooth behavior and the randomness in the model were decoupled using the proposed POD-PCE surrogate model. The time-dependent non-smooth behavior was modeled using the POD modes and the randomness in the model was propagated using the PCE model.

A nonlinear impact oscillator has been investigated for UQ with two different conditions: single contact between the structure and the projectile, multiple contacts between the structure and the projectile. UQ was performed for projectile displacement, projectile velocity and contact force. The stochastic responses were computed very well by the PCE and the POD-PCE models. Effectively, the accuracy of the predicted mean was much high as compared to the predicted standard deviation by both the surrogate models. However, non-physical contact forces were predicted for the multiple impact oscillator by both the surrogate models at some sample points due to the numerical issue with the PCE model which is the main limitation of the proposed approach. For the PCE model, the PCE coefficient vector was identified at each time-step, whereas the stochastic responses were computed using much less POD modes by the POD-PCE model which leads to computing significantly less number of PCE coefficient vector. Hence, the number of PCE model identification was much less for the POD-PCE model without losing the accuracy. As a result, the computational cost for the POD-PCE model was much lower than the other approaches. For a high-dimensional problem, obtaining the POD modes may be inefficient by the POD approach. In that case, the Krylov subspace [START_REF] Freund | Krylov-subspace methods for reduced-order modeling in circuit simulation[END_REF] would be a good alternative to approximate the reduced modes efficiently. Further, the number of PCE coefficients could be reduced by using the sparse PCE model [START_REF] Blatman | Adaptive sparse polynomial chaos expansion based on least angle regression[END_REF] in conjunction with the POD approach.

Fig. 1 :

 1 Fig. 1: 2-DOF nonlinear impact dynamic oscillator

Fig. 3 : 8 Bhattacharyya

 38 Fig.3: Time-dependent statistical moments of the stochastic response quantities for the single impact oscillator: MCS was performed using N MCS = 10 4 model evaluations

Fig. 5 :

 5 Fig. 5: Comparison of responses for single impact oscillator at 5 samples: m p ∈ {342.71, 292.84, 365.62, 310.35, 373.82}g, k c ∈ {13.97, 13.87, 14.61, 14.23, 13.63} × 10 3 MN m -3/2 , v 0 ∈ -{11.32, 9.51, 9.96, 9.11, 10.95}m s -1

Fig. 6 :

 6 Fig.6: Evaluation of mean relative error for all the stochastic response quantities of the multiple impact oscillator

  (N=50,P=2) POD-PCE (N=50,P=3) (a) Mean of f c (zoom on x axis) Standard deviation of f c (zoom on x axis) (N=50,P=2) POD-PCE (N=50,P=3) (f) Standard deviation of ẋp (zoom on x axis)

( a )Fig. 8 :

 a8 Fig. 8: Time-dependent errors for the statistical moments of the response quantities of multiple impact oscillator

Fig. 9 :

 9 Fig. 9:Comparison of responses for the multiple impact oscillator at 5 samples:m p ∈ {342.71, 292.84, 365.62, 310.35, 373.82}g, k c ∈ {13.97, 13.87, 14.61, 14.23, 13.63} × 10 3 MN m -3/2 , v 0 ∈ -{11.32, 9.51, 9.96, 9.11, 10.95}m s -1

Table 1 :

 1 Parameters of the input random variables for the impact dynamic oscillator

	Random variable	Lower bound	Upper bound	Unit
	m p	330 -33	√ 3	330 + 33 √ 3	g
	k

c 16 × 10 3 -1600 √ 3 16 × 10 3

Table 2 :

 2 Accuracy of different surrogate models in assessing the stochastic response quantities for the single impact oscillator by N = 20

				Contact force	Projectile displacement		Projectile velocity
	Method	P n b	n tot	ε	n b	n tot	ε	n b	n tot	ε
	PCE	2	-10010 6.21 × 10 -5	-10010 5.41 × 10 -6	-10010 5.34 × 10 -6
	PCE	3	-20020 14.90 × 10 -5	-20020 2.80 × 10 -6	-20020 7.52 × 10 -6
	POD-PCE 2	5	50	5.54 × 10 -5	2	20	5.53 × 10 -6 3	30	19.40 × 10 -6
	POD-PCE 3	5	100	5.28 × 10 -5	2	40	2.92 × 10 -6 3	60	18.50 × 10 -6

Table 3 :

 3 Efficiency of different approaches for the single impact oscillator

				Contact force Projectile displacement Projectile velocity
	Method	N	P		CPU time (s)
	PCE	20 3	0.061	0.061	0.060
	POD-PCE 20 3	0.003	0.002	0.002
	MCS	10 4 -	41.510	41.510	41.510
	are shown in Figure

  . It seems that both the approaches have predicted

								10 -3
	0.25 0.3				PCE (P=2) PCE (P=3) POD-PCE (P=2) POD-PCE (P=3)			2.8 2.6				PCE (P=2) PCE (P=3) POD-PCE (P=2) POD-PCE (P=3)
								2.4
	0.2							2.2
								2
	0.15							1.8
								1.6
	0.1							1.4
	35	40	45	50	55	60		30	35	40	45	50	55	60
		Number of sample points (N)						Number of sample points (N)
		(a) ε for contact force					(b) ε for projectile displacement
				30	35	40	45	50	55	60
					Number of sample points (N)
							11					Bhattacharyya, VIB-19-1536

Table 4 :

 4 Accuracy of different surrogate models in assessing the stochastic responses for the multiple impact oscillator using

	N = 50
				Contact force	Projectile displacement		Projectile velocity
	Method	P n b	n tot	ε	n b	n tot	ε	n b	n tot	ε
	PCE	2	-	30010 1.14 × 10 -1	-30010 2.52 × 10 -3	-30010 3.91 × 10 -3
	PCE	3	-	60020 1.05 × 10 -1	-60020 1.49 × 10 -3	-60020 2.21 × 10 -3
	POD-PCE 2 31	310	1.13 × 10 -1 3	30	2.52 × 10 -3 9	90	3.91 × 10 -3
	POD-PCE 3 31	620	1.05 × 10 -1 3	60	1.49 × 10 -3 9	180	2.16 × 10 -3

Table 4 .

 4 It is noticeable that almost the same overall accuracies have been achieved by both the surrogates. However, for the POD-PCE model, the PCE coefficient vector was computed only 9 times: it is 3001 for the PCE model.

	13	Bhattacharyya, VIB-19-1536

Table 5 :

 5 Efficiency of different approaches for the multiple impact oscillator Contact force Projectile displacement Projectile velocity

	Method	N	P		CPU time (s)
	PCE	50 3	0.181	0.180	0.181
	POD-PCE 50 3	0.014	0.010	0.014
	MCS	10 4 -	104.950	104.950	104.950
	6 Conclusion

Bhattacharyya, VIB-19-1536