We propose a design methodology to exploit adaptive nanodevices (memristors), virtually immune to their variability. Memristors are used as synapses in a spiking neural network performing unsupervised learning. The memristors learn through an adaptation of spike timing dependent plasticity. Neurons' threshold is adjusted following a homeostasis-type rule. System level simulations on a textbook case show that performance can compare with traditional supervised networks of similar complexity. They also show the system can retain functionality with extreme variations of various memristors' parameters, thanks to the robustness of the scheme, its unsupervised nature, and the power of homeostasis. Additionally the network can adjust to stimuli presented with different coding schemes.

I. INTRODUCTION

anotechnology has long been announced as a game changer of microelectronics. Nanodevices are compact, low power, and above all they provide novel functions like memristivity [START_REF] Strukov | The missing memristor found[END_REF]. However, introducing them into actual products has proved challenging in particular due to their less desirable features: they have strong variability, and it is difficult to achieve acceptable yield. Fully exploiting nanoelectronics potential will thus require architectures and design approaches with an actual immunity to device variability and a high level of fault tolerance.

Spiking neural networks can provide a serious lead since the brain -in particular -itself relies on variable and unpredictable neurons and synapses [START_REF] Marder | Variability, compensation and homeostasis in neuron and network function[END_REF] and manages computational efficiency that outperforms manmade systems. Architectures allying memristors synapses and CMOS neurons have been proposed [START_REF] Snider | Instar and outstar learning with memristive nanodevices[END_REF]- [START_REF] Snider | From Synapses to Circuitry: Using Memristive Memory to Explore the Electronic Brain[END_REF], promising cognitive computing and robustness to defects. In this paper system simulations introduce quantitative results in terms of computing performances and robustness to variability. We show that a simplified spike timing dependent plasticity scheme adapted to memristors is key to achieve effective learning with extreme robustness to memristors' variability. It is associated with the use of unsupervised learning, and of a homeostatic-type mechanism.

We describe the required technology and architecture (section II) and perform system-level simulations on the Manuscript received February 10, 2011. This work was supported by the European Union through the FP7 Project NABAB (Contract FP7-216777).

D. Querlioz is with Institut d'Electronique Fondamentale, Univ. Paris-Sud, CNRS, 91405 Orsay, France (phone: +33 1 61 15 33 58; fax: +33 1 69 15 40 20; e-mail: damien.querlioz@u-psud.fr).

O. Bichler and C. Gamrat are with CEA, LIST, Embedded Computers Laboratory, 91191 Gif-sur-Yvette Cedex, France (e-mail: olivier.bichler@cea.fr, christian.gamrat@cea.fr)

MNIST database [START_REF] Lecun | Gradient-based learning applied to document recognition[END_REF] that show the potential of the approach and its robustness. The system consists of a first fully unsupervised layer that extracts features of the inputs using a simplified spike timing dependent plasticity (section III), and can be complemented with a layer to label the outputs (section IV). The network performance compares favorably with traditional -but supervised networks with similar numbers of adjustable parameters and achieves excellent tolerance to various memristors' parameters variability.

Technological context: In recent years several classes of adaptive devices have emerged. The most famous are memristors [START_REF] Strukov | The missing memristor found[END_REF] (stricto sensu memristive devices), resistive RAMs [START_REF] Burr | Overview of candidate device technologies for storage-class memory[END_REF], or adaptive transistors like the NOMFETs [START_REF] Alibart | An Organic Nanoparticle Transistor Behaving as a Biological Spiking Synapse[END_REF] and OG-CNTFETs [START_REF] Agnus | Two-Terminal Carbon Nanotube Programmable Devices for Adaptive Architectures[END_REF]. An emerging route is to exploit them as synapses in neuromorphic circuits (electronic circuits that work analogously to the brain). In particular it has been suggested [START_REF] Snider | Spike-timing-dependent learning in memristive nanodevices[END_REF], [START_REF] Pérez-Carrasco | On neuromorphic spiking architectures for asynchronous STDP memristive systems[END_REF], [START_REF] Afifi | Implementation of biologically plausible spiking neural network models on the memristor crossbar-based CMOS/nano circuits[END_REF], and shown experimentally [START_REF] Jo | Nanoscale Memristor Device as Synapse in Neuromorphic Systems[END_REF], that such devices could indeed reproduce a feature of biological synapses -spike timing dependent plasticity (STDP) [START_REF] Markram | Regulation of Synaptic Efficacy by Coincidence of Postsynaptic APs and EPSPs[END_REF], [START_REF] Bi | Synaptic modification by correlated activity: Hebb's Postulate Revisited[END_REF] -that is believed to be a foundation for learning in the brain [START_REF] Dan | Spike Timing-Dependent Plasticity of Neural Circuits[END_REF], [START_REF] Gorchetchnikov | A model of STDP based on spatially and temporally local information: Derivation and combination with gated decay[END_REF]. A system consisting of nanoscale synapses and CMOS processing units (neurons) could be a major breakthrough in computing, possibly allowing cognitive-type tasks. This precise idea is currently receiving considerable interest through for example the DARPA SyNAPSE program in the USA or similar projects in Europe and Asia. However, its sustainability is still to be demonstrated.

As mentioned above, one common characteristic of all these technologies is a high degree of variability [START_REF] Jo | CMOS Compatible Nanoscale Nonvolatile Resistance Switching Memory[END_REF], [START_REF] Borghetti | Memristive' switches enable `stateful' logic operations via material implication[END_REF]. Though technological maturity will improve the situation, it is expected to be intrinsic especially if we scale to ultra-small sizes that make the devices particularly sensitive to any variation in the fabrication process. This puts a challenge to designing systems able to exploit their new functionality. This paper provides some insight on how such technology may be exploited.

Prior art: Other proposals have been made to exploit variable adaptive devices in the context of nanotechnological implementations. Most proposed architectures rely on reconfigurable logic (FPGA type) [START_REF] Heath | A defecttolerant computer architecture: Opportunities for nanotechnology[END_REF], [START_REF] Dehon | Array-based architecture for FET-based, nanoscale electronics[END_REF] or on state-based supervised neural networks [START_REF] Lee | Defect-tolerant nanoelectronic pattern classifiers[END_REF]. In the first approach variability is addressed through error mapping and redundancy, in the second case through a standard neural network approach using supervised learning based on error gradient descent.

Our approach of using unsupervised learning with asynchronous spiking neural networks to tackle the

Simulation of a Memristor-Based Spiking Neural Network Immune to Device Variations

Damien Querlioz, Olivier Bichler and Christian Gamrat N variability issue of nanodevices is original and takes inspiration from recent works in computational neuroscience and neural networks [START_REF] Guyonneau | Neurons Tune to the Earliest Spikes Through STDP[END_REF]- [START_REF] Nessler | STDP enables spiking neurons to detect hidden causes of their inputs[END_REF]. Different works have been published that go into that direction. As mentioned above, several proposals exist to use memristive devices for STDP [START_REF] Snider | Spike-timing-dependent learning in memristive nanodevices[END_REF], [START_REF] Pérez-Carrasco | On neuromorphic spiking architectures for asynchronous STDP memristive systems[END_REF], [START_REF] Afifi | Implementation of biologically plausible spiking neural network models on the memristor crossbar-based CMOS/nano circuits[END_REF]. In this paper we propose a simplified STDP scheme that we believe to be easier to implement, and to allow learning in a more natural way. Additionally, we propose a homeostasic property to the neurons that has not been used in this context before and that is shown to be essential for the robustness of the scheme. One work had already shown that memristors with STDP could allow the emergence of receptive fields in a variability-compatible unsupervised approach and synchronous neurons [START_REF] Snider | Self-organized computation with unreliable, memristive nanodevices[END_REF]. Our work uses asynchronous designs (like the ones used in the neuromorphic community [START_REF] Arthur | Learning in silicon: Timing is everything[END_REF], [START_REF] Livi | A current-mode conductance-based silicon neuron for address-event neuromorphic systems[END_REF]), and performs full learning on a standard dataset. Finally, an alternative way to allow learning with memristors in a variability-compatible way can be to use all digital designs [START_REF] Snider | From Synapses to Circuitry: Using Memristive Memory to Explore the Electronic Brain[END_REF]. This requires more devices per synapses [START_REF] Snider | Instar and outstar learning with memristive nanodevices[END_REF]. In this paper we show that variation-tolerance can be retained by using nanodevices with continuous variation of the conductance.

II. PRESENTATION OF THE NETWORK AND OF ITS IMPLEMENTATION

A. Architecture

In this paper we first propose a simple two layer topology. Input and output CMOS neurons are connected in a feedforward manner by nanodevices (the synapses), using a crossbar configuration (described below). The synapses learn using a simplified STDP scheme (II.B.1). The output neurons behave as leaky integrate-and-fire neurons (II.B.2) and have a homeostatic property (II.B.4). They are also connected by inhibitory connections using diffuser networks (II.B.3)

The input neurons present the stimuli as asynchronous voltage spikes using several possible coding schemes described in section III (spiking rate is proportional to stimulus intensity). These stimuli may originate for example directly from a spiking retina [START_REF] Lichtsteiner | A 128× 128 120 dB 15 μs Latency Asynchronous Temporal Contrast Vision Sensor[END_REF] or cochlea [START_REF] Chan | AER EAR: A Matched Silicon Cochlea Pair With Address Event Representation Interface[END_REF] designed in the neuromorphic community. It is natural to lay out the nanodevices in the widely studied crossbar as illustrated on Fig. 1 (CMOS silicon neurons and their associated synaptic driving circuitry are the dots, the squares being the nanodevices). The synapses indeed act as adaptive resistors. With the crossbar layout, if several synapses are active at the same time (i.e. receiving spikes), the output receives directly the sum of the currents flowing through the synapses. In a more futuristic design, the system could also be laid out in a CMOL architecture where nanodevices crossbar is fabricated on top of the CMOS neurons and driving circuits [START_REF] Lee | Defect-tolerant nanoelectronic pattern classifiers[END_REF].

As a result of learning, the output neurons will become selective to the different stimuli classes that are presented in a fully unsupervised manner: the output neurons will develop selectivity to specific features contained in the input patterns. The learning rule of the nanodevices needs to be fully local to be implementable in the crossbar architecture. The behavior of the neurons needs to be simple to make it easy to implement in a compact way. We now describe how this can be achieved.

B. Neurons and synapses

1) Synaptic behavior: In this scheme, the synapses are acting in two ways. They are variable resistors and thus transmit spikes with a variable conductance (or weight). But Fig. 2. Pulses for simplified STDP (voltage pulses as a function of time). When an input neuron spikes, it applies a PRE pulse to the nanodevices to which it is connected. When an output neuron spikes it applies a POST spike to the nanodevices to which it is connected. When the voltage applied on the device (difference between the voltages applied at the two ends) reaches VT+ or VT-, its conductance is increased or decreased, respectively. they also adapt their weight depending on the activity of the neurons to which they are connected, which provides the foundation of learning by the system.

The memristive nanodevices are programmed as follows. When applied a positive voltage pulse higher than a threshold V T+ , they increase their conductance. When applied a voltage pulse smaller than a negative threshold V T-, they decrease their conductance [START_REF] Jo | Nanoscale Memristor Device as Synapse in Neuromorphic Systems[END_REF], [START_REF] Snider | Spike-timing-dependent learning in memristive nanodevices[END_REF]. Previous works have shown that memristors can implement spike timing dependent plasticity (STDP), a learning rule used by synapses in the brain. This was proposed in [START_REF] Pérez-Carrasco | On neuromorphic spiking architectures for asynchronous STDP memristive systems[END_REF], itself close to the proposal [START_REF] Snider | Spike-timing-dependent learning in memristive nanodevices[END_REF] and experimentally demonstrated in [START_REF] Jo | Nanoscale Memristor Device as Synapse in Neuromorphic Systems[END_REF]. We propose a simpler adaptation of these schemes (Fig. 2).

• When an input neuron spikes it applies a long pulse to its synapse (PRE pulse). This voltage is high enough to drive some current into the memristor but not enough to reprogram it. This current is integrated by the output neurons (resistor role of the synapse). If several synapses connected to the same output neuron are ON at the same time their currents are summed. • When an output neuron spikes it applies a pulse that is a succession of a negative and a positive bias (POST pulse). If no PRE pulse is being applied on the device, only the second part reaches a threshold and the weight of the synapse is decreased by w δ -. However, if the input neuron had spiked recently, and the PRE pulse is still being applied on the other end of the device, the voltage applied on the device actually increases its weight by w δ + , as illustrated on Fig. 2. This simple learning rule, easily implemented with nanodevices, is the ground of all learning. Compared with the purely bioinspired and more complex scheme introduced in [START_REF] Pérez-Carrasco | On neuromorphic spiking architectures for asynchronous STDP memristive systems[END_REF], no delay matching is necessary between the PRE and POST synaptic waveforms, which should make the driving circuitry much easier to design. The learning rule is also more appropriate for our feed-forward application.

The way this simple learning rule works is straightforward. When an output neuron declares a spike (at time t spike) , it increases by w δ + the weights of the synapses connected to input neurons that spiked recently (from t spike to t spike -t PRE , if t PRE is the duration of the PRE pulse). It decreases by w δthe weights of all its other synapses. This increases the sensitivity of the neuron to the specific pattern that activated it, making it more likely to spike for a similar (correlated) pattern in the future. This process -that works surprisingly well in practice, as we show in this paper -has been partially theorized in [START_REF] Nessler | STDP enables spiking neurons to detect hidden causes of their inputs[END_REF].

To model the weight increments and decrements w δ + and w δin our system simulations, we take inspiration from the memristor measurements [START_REF] Jo | Nanoscale Memristor Device as Synapse in Neuromorphic Systems[END_REF], [START_REF] Pickett | Switching dynamics in titanium dioxide memristive devices[END_REF]. When a POST occurs: if PRE has occurred just before the weight is increased by

The exponential factor expresses the fact that is observed in most memrsitive technology that a same voltage pulse has reduced effect on the device conductance if applied several times [START_REF] Jo | Nanoscale Memristor Device as Synapse in Neuromorphic Systems[END_REF], [START_REF] Pickett | Switching dynamics in titanium dioxide memristive devices[END_REF]. This natural (partly multiplicative) device behavior is actually a valuable behavior for learning. The parameters α + , α -, β + , β + depend heavily on the PRE and POST pulse voltages that are chosen. These parameters, as well as minimum and maximum weights w min and w max are subject to device variability.

2) Output neurons' dynamics: Exploiting the devices requires connecting them to processing units -silicon neurons able to process and generate spikes in a bioinspired manner by integration of their input. An approach widely studied in the neuromorphic community is to use analog circuits (generally with transistors operating in the subthreshold regime) able to receive and generate asynchronous spikes [START_REF] Arthur | Learning in silicon: Timing is everything[END_REF], [START_REF] Livi | A current-mode conductance-based silicon neuron for address-event neuromorphic systems[END_REF]. Though smaller than in nanodevices, variability is also a serious problem for such devices. It is a challenge for any neuromorphic design [START_REF] Chan | AER EAR: A Matched Silicon Cochlea Pair With Address Event Representation Interface[END_REF] and will become even more crucial when scaling to modern technology.

Neurons are leaky integrate-and-fire type that is meant to solve the simple equation (expressed in normalized unit):

input I gV dt dV = + τ (3)
where V is the state variable (a current or a voltage) of the neuron. I input is the current flowing through the line of the crossbar connected to the neuron:

∑ = j j input I I ( 4
)
where the I j are the currents flowing through each memristor j connected to output neuron.

The neuron declares a spike if V reaches a given threshold V th , in which case V is reset to zero and the integration is prevented during a so-called refractory period.

Such neurons have been implemented with compact CMOS neuromorphic circuits as in [START_REF] Indiveri | A VLSI array of lowpower spiking neurons and bistable synapses with spike-timing dependent plasticity[END_REF] or equivalently with a neuromorphic circuit with positive feedback as in [START_REF] Arthur | Learning in silicon: Timing is everything[END_REF], [START_REF] Livi | A current-mode conductance-based silicon neuron for address-event neuromorphic systems[END_REF].

This kind of CMOS design is extremely low power because most transistors operate in the subthreshold regime, and thanks to the use of asynchronous computation [START_REF] Arthur | Learning in silicon: Timing is everything[END_REF], [START_REF] Livi | A current-mode conductance-based silicon neuron for address-event neuromorphic systems[END_REF]. On the other hand, power consumption of the nanosynapses will depend tremendously on the nanotechnology, and should scale with the devices.

3) Output neurons' lateral inhibition: When an output neuron spikes it should send inhibitory signals to the other output neurons of the layer that prevent them from spiking during the inhibition time and resets their potential to zero.

This inhibition between the neurons can be implemented in a compact way through diffuser networks as in [START_REF] Arthur | Learning in silicon: Timing is everything[END_REF]. With this inhibition, the network is reminiscent of a Winner-Takes-All topology [START_REF] Nessler | STDP enables spiking neurons to detect hidden causes of their inputs[END_REF].

More precisely, when an output neuron spikes, the state variable V of the other output neurons is reset to zero during a time t inhibit .

4) Homeostasis: A final difficulty for the architecture is the adjustment of the neurons' threshold. It is essential, but although this is natural for traditional artificial neuron networks, it is not for spiking neurons. A bioinspired original route is homoeostasis [START_REF] Marder | Variability, compensation and homeostasis in neuron and network function[END_REF]. A target activity is defined for the neurons (i.e. a number of times an output neuron should spike over an extended period of time, like 100 digits presentation). Regularly the threshold of the neuron is increased if the average activity of the neuron is above the target, and decreased if it is below.

This ensures that all the outputs neurons are used and adjust the neurons' thresholds to the stimuli for which they become specialized. In neuromorphic hardware, homoeostasis has been implemented with analog memories like in [START_REF] Bartolozzi | Implementing homeostatic plasticity in VLSI networks of spiking neurons[END_REF] or could be implemented digitally. The potential associated with this technology is evidenced in this paper in section III.B.

C. Simulations

In this paper, all simulations are system-level and are based on a C++ special purpose code. The code is eventbased for simulation performance and runs on traditional CPUs. Simulation parameters as introduced above are τ = 100 ms, g = 1, V th = 0.5 (normalized unit), inhibition time t inhibit = 10 ms, α + = 10 -2 , α -= 5.10 -3 , min w = 10 -4 , max w = 1 (normalized units), β + = 3.0, β -= 3.0. The width of the PRE pulses is 25 ms. Parameter variations are introduced around all the parameters using Gaussian random numbers (the value of their standard deviation is given in Section III). The initial weights are selected randomly around mid-range (0.5). The stimuli are applied using the coding schemes described in section III.C.

For demonstration of the concept, in this paper we use the widely studied case of handwritten number recognition. The MNIST database is used, which consists in handwritten number of 28x28 pixels by 250 writers [START_REF] Lecun | Gradient-based learning applied to document recognition[END_REF]. This database has been used as test for many learning algorithms.

In order to achieve learning, we present the full MNIST training database (60,000 digits) three times to the circuit. Each input neuron is connected with one pixel of the image. It emits spikes with a jittered rate that is proportional to the pixel intensity (maximum rate is 20 Hz). The initial phase is random. Every digit is presented during 350ms. No kind of preprocessing on the digits is used and the set is not augmented with distortions. The network is then tested on the MNIST test database (10,000 digits).

Fig. 3 plots the synaptic weights learned by the system in a configuration with only 10 output neurons. It is remarkable that without any supervision and using only our local custom STDP rule, the system has identified 9 (out of 10) different numbers, the real features of the input. Moreover it has learned the distinctive features of the digits (and not just the most likely handwriting): it has learnt the loop of the digit two, the bottom of the six, or the horizontal parts of three and eight.

In order to evaluate quantitatively the network performance, Fig. 4 plots the recognition rate during the learning process for different numbers of output neurons. To evaluate the recognition rate, we associate output neurons with the digit for which they spike the most frequently a posteriori. In hardware this association could be performed with digital circuitry. This labeling could also be implemented using nanodevices in a variability compatible way as will be shown in section IV.

With ten output neurons the recognition rate reaches 60 %. With 50 output neurons it reaches 81% and with 300 output neurons 93.5 % (obtained with the same number of adjustable weights). A traditional artificial neural network with back-propagation with 300 hidden neurons reaches 95 % [START_REF] Lecun | Gradient-based learning applied to document recognition[END_REF], which compares to our rate of 93.5 %. In the literature, the best algorithm has a largely superior 99.7 % recognition rate, but has 12 million adjustable weights (vs. 235,200 here) and is using a largely augmented training set [START_REF] Ciresan | Deep Big Simple Neural Nets Excel on Handwritten Digit Recognition[END_REF]. Though our numbers are clearly more modest, the interest here is that the network is fully unsupervised with local extremely simple learning rules and variability immunity.

III. VARIABILITY IMMUNITY

A. Synaptic variability

We first study the impact of nanodevices (synaptic) variability on the performance of the network. For this study, we simulated the network with 50 output neurons, introducing different kinds of variability. All the simulations were repeated five times, and the recognition rate given is averaged on the five runs. Simulation time was about 8 hours per run on an AMD Opteron 2216 CPU. Results are reported in Table I.

On the first row of Table I, we evaluate the impact of variability of the initial state (i.e. the initial weight of the nanodevices, before learning starts). We can see it has no impact on the final performance, even when the variability is extreme. This means that we do not need to control the initial weight precisely. This is fortunate since this would be extremely difficult with memristors.

On the second row of Table I, we evaluate the impact of the variations on the learning increments and decrements ( α + , α -). These variations can be caused by variations of the device thresholds and of their programming characteristics. We can see that the network is immune to variations up to 25% of the mean value of these parameters, which is already an extremely high level of variation for an electron device (but typical for research nanodevices). With 50% of variation, there is a small decrease of performance. With a dramatic variability of 100% on the synaptic parameters, the recognition rate decreases significantly, but interestingly the functionality of the network is not challenged. We should note that with a variability of 50%, 4% of the nanodevices cannot be programmed in at least one direction (i.e. they have an α value of 0). The latter figure becomes 30% when parameter variability is 100%. This is an indication of the overall tolerance to defects.

The last row of Table I adds variations of min w , max w . In that case some nanodevices will overpower the others due to increased maximum weight. Results are however similar to the previous case, and extreme robustness is preserved. The most sensitive parameter is the maximum weight.

Tolerance to synaptic variability is natural in a way since learning actually takes its root in synaptic variability. The degree of tolerance achieved is however still impressive.

B. CMOS neuron variability

We now study the impact of the variability of the CMOS neurons. The impact of the variability of their threshold is seen in Table II. We can see that without homeostasis the impact of the variability is dramatic. The more excitable neurons (the ones with lower thresholds) spike predominantly and the other neurons do not specialize efficiently. In a typical simulation with a threshold variation of 25 %, the most excitable neuron spikes 51 % of the time, and in a simulation with a threshold variation of 50 %, it spikes 91% of the time. Homeostasis however fully compensates this issue (all the neurons spike between 1.5 and 3% of the time, most neurons spiking around 2%). And the same degree of performance than without threshold variation is achieved as evidenced in Table II. Homeostasis appears as extremely valuable in this application.

C. Stimulus encoding

Another remarkable point is the insensitivity to stimulus encoding scheme. Different schemes are possible to encode the stimulus into spikes as illustrated in Fig. 5. In the first two schemes, the input neurons spike periodically, the firing rate being proportional to the input quantity. The input neurons can spike either in phase or out of phase. This is a natural scheme based on how spiking retinas currently work [START_REF] Lichtsteiner | A 128× 128 120 dB 15 μs Latency Asynchronous Temporal Contrast Vision Sensor[END_REF]. In the last possibility, input neurons spike in a Poissonian way, the time constant being inversely proportional to the input. All three coding styles lead to a performance between 81% and 82%, which again shows the robustness of the approach. This kind of robustness to spike jitter is consistent with computational neuroscience works on STDP [START_REF] Guyonneau | Neurons Tune to the Earliest Spikes Through STDP[END_REF].

IV. OVERALL NETWORK

A. Supervised learning scheme

The second layer of the proposed architecture is a simple classical supervised classifier. It is used to label outputs of the first unsupervised layer.

For this second layer, we propose a new scheme. It is similar to the unsupervised scheme presented in the previous sections with the difference that no voltage is applied on the nanodevices when an output neuron spikes. Instead, after a pattern has been applied for a given time, input and output neurons apply a programming pulse on the nanodevices as illustrated on Fig. 6.

The input neurons apply a pulse only if they spiked during the period. Depending if the output neuron reacted as was expected, it will program the nanodevices connected to input neurons that spiked one way or another. The pulse from the input neurons (PRE) indeed contains a positive and a negative part so that by conjunction with the pulse applied to the output neuron (POST) it may program the nanodevice one way or another.

If the output neuron was expected to spike (it corresponds to the digit actually presented), the connections to the input neurons that spiked should be strengthened. The output neuron applies a negative spike that will strengthen the positive part of pre-synaptic pulses, and increase the weights of nanodevices that present such pulses. If the output neuron was not expected to spike, inversely it should apply a positive pulse to weaken the synapses connected to input neurons that spiked. If the output neuron actually spiked (i.e. made a mistake) a significant pulse is applied to significantly reduce the weights. If it did not spike (and thus was right)which is the most typical situation -the pulse is smaller in order to only weakly weaken the weights (and not cause the weights to decrease too rapidly).

B. Two layer scheme

The role of the supervised layer is to label neurons trained by the unsupervised first layer of section III. As illustrated on Fig. 7, this works well and in that case the insensitivity to variations is kept because labeling is not a difficult problem. Even with variations of α + , α -, min w , max w of 50% the supervised layer performance is matching perfectly the unsupervised layer.

The supervision does not improve performance, but provides a practical approach for the labeling operation. This opens up many opportunities. One promising way could be to pretrain the unlabeled layer on natural data, and then to perform a "labeling" on a small training set. This could be a useful way to process all kinds of natural data.

V. CONCLUSION

In this paper we have shown how by using a simple custom STDP scheme, memristors associated with CMOS neuromorphic circuits could perform unsupervised learning in a way that is extremely robust to variability. Nonsupervision is the foundation for this tolerance and will provide extreme adaptability to the system. Homeostasis was also introduced as an essential component to retain insensitivity to CMOS variability. In the last section a scheme to use memristors for supervised learning was also introduced in a second layer. When used for labeling results of the first unsupervised layer, it retains excellent variability robustness, and provides a scheme for complete learning of real case data. This design approach could be the groundwork for future circuits that will process natural data in a compact and low power way. Thanks to their unsupervised learning capabilities, such circuits will be able to adapt to various environments. Future work should focus on the experimental demonstration of these concepts, and to demonstrate its scaling to more complex multi-layer networks, and to other kinds of sensory stimuli.

Fig. 1 .

 1 Fig. 1. Circuit topology. Wires originate from CMOS input layer (horizontal black wires) and from the CMOS output layer (vertical grey wires). Adaptive nanodevices are at the intersection of the horizontal and vertical wires.

Fig. 4 .Fig. 3 .

 43 Fig. 4. Recognition rate during learning for simulations with different numbers of output neurons (from bottom to top: 10, 50, 100, 300). The recognition rates are running averages on 10,000 iterations.

Fig. 7 .Fig. 6 .

 76 Fig. 7. Two layer learning. Learning rate of the second (labeling) layer trained by variability compatible supervision (red) compared with learning of the first layer (black), as computed in the previous section. Variability of 25% on all parameters. (initial weights, α + , α -, min w , max w )

Fig. 5 .

 5 Fig.5. Illustration of the three stimulus coding schemes. Stimulus spikes timing for ten input neurons with the same pixel value in the three schemes. All three lead to a recognition rate between 81% and 82%, without any parameter needing to be adjusted.

TABLE I RECOGNITION

 I RATE ON THE MNIST TESTING SET WITH DIFFERENT KINDS OF NANODEVICE (SYNAPTIC) VARIABILITY. ALL THE SIMULATIONS WERE REPEATED 5 TIMES, THE RECOGNITION RATE GIVEN IS AVERAGED ON THE 5 RUNS. ALL THE SIMULATIONS HAVE HOMEOSTASIS. THE DEGREE OF VARIABILITY (10%, 25%, 50%, 100%) IS DEFINED AS STANDARD VARIATION OF THE MEAN VALUE.

	Dispersion on the parameter	10%	25%	50%	100%
	Initial weights	81.3% 82.0% 80.8% 81.3%
	Learning steps α + ,	81.8% 81.4% 79.0% 74.0%
	α -
	Learning steps α + ,	81.9% 80.6% 77.2% 67.8%
	α -, min. and max.
	weights w min , w max

TABLE II IMPACT

 II OF THE NEURONS' THRESHOLD VARIABILITY, WITH AND WITHOUT HOMEOSTASIS ON THE RECOGNITION RATE. ALL THE SIMULATIONS WERE REPEATED 5 TIMES. THE DEGREE OF VARIABILITY (10%, 25%, 50%) IS DEFINED AS STANDARD VARIATION OF THE MEAN VALUE.

	Variability on the neurons' threshold	Homeostasis	No homeostasis
	No variability	81.3%	73.4%
	10%	81.9%	60.2%
	25%	81.2%	27.5%
	50%	80.8%	14.9%

ACKNOWLEDGMENT

The authors would like to thank S. J. Thorpe, W. Zhao, J.O. Klein, J. Saint-Martin, P. Dollfus, D. Deleruyelle and B. De Salvo for fruitful discussions and collaboration.


F Jacquemin

\a Vautrin

R Guilh3n

Internal stress reduction in filament wound composite pipes by humid conditioning : numerical simulation

Keywords: Internal stress, hygrothermal treatment, filament winding process, polymer composite pipes

The paper introduces a new approach which tends to lower the internal stress state within composite pipe walls by using special post-cure hygrothermal treatment. Residual curing stresses in polymer composites are particularly critical in dry conditions, it is proved that an easy way to limit their effect on the material integrity is to introduce hygrothermal stresses. The residual stresses define the initial state of the pipe before it's service life. Usual ambient conditions generate internal stresses which add to the residual manufacturing stresses. It is shown that particularly hygrothennal conditions introduced just after the manufacturing can induce transient stresses which finally will reduce the global stresses. The proposed approach makes it possible to check the influence of timevariable process temperature and pressure, and hygrothermal in-service conditions on polymer matrix composites behavior.

Introduction

The composite pipes are usually manufactured by filament winding process. The filament winding process is commonly composed of three steps : (1) winding the composite tows onto a mandrel, (2) curing, and (3) cool-down and removal of the mandrel. After winding, the assembly composed of the mandrel and the composite pipe is introduced in an autoclave to cure the epoxy resin according to a prescribed temperature and pressure history [ 1]. The stress free state is supposed to be achieved at the cure temperature then there is no stress development prior to cool-down. In thin composites a uniform through-the-thickness temperature distribution assumption is justified. Such an approach is not appropriate for predicting the process-induced stresses in thick thermoset composites where complex temperature develops during the manufacturing process (2]. The residual stresses in the pipe wall, due to non-uniform temperature fields, are computed within the framework of thermoelasticity. Special attention is then paid to the particular effect of the mandrel thickness on the residual stress development. When curing is completed, the mandrel is removed and the pipe is submitted to the environmental conditions. Ambient conditions generate internal stresses which add to the residual manufacturing stresses. Such residual stress growth may damage pipes manufactured by using thick mandrels. The scope of the paper is to show that the introduction of particular hygrothermal conditions just after the manufacturing process induce transient stresses which counteract the residual curing stresses.

Manufacturing process

Temperature field.

We consider a pipe made up of n perfectly bonded plies, whose inner and outer radii are a and b respectively, manufactured by filament winding process onto a mandrel. Any ply i of it is a cylinder whose inner and outer radii are r; and ri+J respectively. The pipe and the mandrel are initially at the same uniform temperature T 0 . Thermal transfers take place on the external face that are controlled by a coefficient of exchange h kept constant. The autoclave temperature T(t) follow the cure cycle. It is assumed that no heat transfer occurs at the inside of the mandrel [START_REF] Wisnom | Proc. ICCM[END_REF]. Conditions of continuity of the temperature and heat flux are prescribed at ply-interfaces.

The temperature field in ply i (the internal ply being the mandrel) is solution of the equation set bellow, with Fourier's equation (1) and boundary and initial conditions (2) :

2 ori a Tj 1 aTi . - = k • ( --+--) a< r < b 1 = 1 to n at I Or2 r Or ' ' (1)
Ti(rj,t) = Ti+I(rj,t)

A• 8Ti (ri, t) _A• 8Ti+I (ri, t)

I ar -t+I ar 8T 1 (a,t)=O ar (2) aTn~b, t)-h(T(t)-Tn (b, t)) = 0 Ti (r,O) =To
The thermal field, in the composite pipe during the cool-down, is determined by using a finite difference scheme.

Process-induced stresses.

The residual stresses for every ply at any time are calculated by using the classical equations of solid mechanics : constitutive laws of thermoelastic orthotropic materials (3), strain-displacement relationship, compatibility and equilibrium equations and boundary conditions [START_REF] Jacquemin | [END_REF].

cr = L: (E-a(T-To)) (3)
At the beginning of the cool-down, let us consider a stress-free state. Moreover a and L, respectively the tensors of thermal expansion coefficients and stiffnesses, are assumed to be constants.

The axial and circumferential components of the displacement field of the ith ply, respectively u and v, are then express ( 4) :

(4)

The radial component of the displacement field w satisfies the following equation :

The temperature field is determined by using a finite difference scheme. To propose a close form solution of the displacement, we subdivide the mandrel-composi te structure and assume on each subdivision a parabolic temperature field (6) :

2 T =A•r +B•r+C• I I I I (6)
Thus, we obtain the radial component of the displacement field, solution of the equation ( 5), for each subdivision : (7) Finally for n-layer tube the displacement field depends on 4n constants (R 1 , R 2 , R 3 and~) to be determined. The constants are calculated from the following conditions : continuity of the displacement components at each interply; global force balance of the cylinder; continuity of the normal stress at each interply and its nullity on the internal surface and an autoclave pressure applied on the external surface.

Internal stresses generated by the ambient conditions

When curing is completed, the mandrel is removed and the pipe is submitted to the environmental conditions. Ambient conditions generate internal stresses which add to the residual manufacturing stresses. The internal stresses in the pipe wall are computed within the framework of hygro-thermo-ela sticity. We consider an orthotropic hygrothermoelast ic constitutive law, straindisplacement relationship and compatibility and equilibrium equations [START_REF] Jacquemin | [END_REF]. To calculate the unknown constants, 4n constants for n-layer tube, we applied identical relations to the determination of the process-induced stresses but with the nullity of the normal internal stress on the two lateral surfaces of the composite pipe.

Example

We consider a thick walled carbon-epoxy (T300/5208) pipe made up of 20 plies of equal thickness, the inner and outer radii being respectively 75 mm and 90 mm. The lamination sequence (corresponding to high buckling strength) studied is [9061'60/-60/(30/-30) 5 /60/-60/902]. An aluminium mandrel 6061-T6 is used for the manufacturing process [5]. The mechanical properties of the composite (T300/5208) and aluminium mandrel are respectively presented in Tables 1 and2.

The cure temperature and room temperature are 180 °C and 20 oc respectively. The cooling rate is 3°C/min. The thermal properties of the composite (T300/5208) and aluminium mandrel are presented in Table 3. The thickness of the mandrel hm is chosen to be equal to one-third, one-fifth, one-seventh or one-ninth of the total thickness h (mandrel+ composite).

Table 1 :Mechanical properties in the reference frame linked to the fibres (

1 ,2,3) Material Et [Gpa] E2, E3 [GPa] v12, Yt 3 v 23 G12 [GPa] a 1 [K-1 ] a 2 , a 3 [K-1 ]
T300/5208 181 10,3 0,28 0,43 7,17 0,02.10-6 22,5.10-6

Table 2: Mechanical properties for aluminium mandrel (6061-T6)

Material E [GPa] v a [K-1 ]
Aluminium 72.5 0,3 23.6.10-6

Table 3 :Thermal properties of the composite (T300/5208) and aluminium mandrel (6061-T6)

Materials 'A [W.m• 1 .K-1 ] p [kg.m-3 ] c [J.kg• 1 .K-1 ] h [W.m• 2 .k-1 ]
T300/5208 0,7 1590 857 50

Aluminium 180 2700 964

The use of the Tsai strength criterion allows us to estimate the internal stresses versus the ply resistance. The tensile failure stresses in the longitudinal, normal and transverse direction denoted respectively X, Y, Z are reported in Table 4. The corresponding failure stresses for compression are X', Y', Z' and S stands for the plane shear failure stress. These strengths are assumed not to depend on humidity and temperature. R, the strength factor derived from the Tsai strength criterion, stands for the ratio between the ultimate stress and the applied stress. Therefore, the failure occurs when R is less or equal to 1. The cool-down from the cure temperature (180 °C) to the room temperature (20 °C) for a cooling rate of3°C/min takes 3200 s. The thermal fields at the end of the cool-down, for different thickness of the mandrel, are plotted in Fig. l(a). There is a considerable lag between the tube and autoclave temperatures, especially at the inner surface where no heat transfer occurs. Thermal gradients are growing functions of the mandrel thickness. Fig. 1(b) shows that the minimum strength factor is reached at the inner surface (90° layer). Even if the strength factor remains encouraging for the composite plies, we recommend the utilization of thin mandrel which one the strength factor is higher on the inner surface. When curing is completed, the mandrel is removed and the pipe is submitted to the environmental conditions : a temperature of 20 oc and a relative humidity of 50 %. Since the temperature diffusion speed exceeds the moisture diffusion one by several orders, it is reasonable to assume that a uniform temperature of 20 oc is reached in the composite pipe without change of the moisture concentration (dry conditions). This change of temperature induced internal stresses which add to the residual manufacturing stresses. Fig. 2(a) shows the corresponding strength factor. Critical values of the strength factor, less to 1, are reached for the pipe manufactured by using thick mandrels : one-third and one-fifth of the total thickness. To reduce the internal stress state a humid conditioning, temperature of 60 oc and a relative humidity of 50 %, is introduced before the placement at ambient conditions of the composite pipe. We check, in Fig. 2(b ), that the change in temperature, from the temperature at the end of the cool-down (Fig. 1) to a uniform temperature of 60 °C, no involves damage in the composite pipe for all the mandrels considered. The change of the moisture concentration from the dry conditions to a relative humidity of 50 % produce compressive stresses which induced a growth of the strength factor (plotted in Fig. 3(a) for values less than 1 0). Fig. 3(b) depicted the strength factor, for values less than 10, when the composite pipe is placed in the environmental conditions. The strength factor undergoes important gradients and it reaches a minimum value of about 2, 75 on the inner surface. Finally, the material integrity is preserved by introducing transient hygrothermal stresses. 1,5 :

6 1,4 ~ v <!! "&, 1,3 c ~ 1,2 ; 0,9 I 0,8 ' -- 75
radius [nm) radius [mm] -3 h m -.--shm ~7hm -+-9hm -3 h m -.--shm -7 h m -+-9hm (a) (b)

Conclusion

The problem of internal stresses induced by processing and ambient conditions in thick polymer composite pipes is examined. The paper present an approach that allows to reduce the internal stress state and prevent process-induced damage by using a post-cure hygrothermal treatment. The growth of internal stresses, induced by the ambient conditions after the manufacturing process, can lead to damage for thick composite pipes. It is shown that the introduction of a humid conditioning, before the placement in the ambient conditions of the composite structure, involves a reduction of the internal stresses. The influence of the mandrel thickness on the residual stress development is investigated. It is emphasized that curing stresses can be manipulated by choosing tailored mandrel thickness.

Figure 1 :Figure 2 :Figure 3 :

 123 Figure 1 : At the end of the cool-down: (a) thermal fields, (b) strength factors

Table 4 :

 4

		T300/5208 strengths [6]
	Strengths [MPa]	X,X '	Y,Z	Y ',Z'	s
	T300/5208	1500	40	246	68


Tamara Al-Bittar
email: tamara.albittar@gmail.com

Abdul-Hamid Soubra
email: abed.soubra@univ-nantes.fr

Efficient Sparse Polynomial Chaos Expansion Methodology for Computationally-Expensive Deterministic Models

The sparse polynomial chaos expansion (SPCE) methodology is an efficient approach that deals with uncertainties propagation in case of high-dimensional problems (i.e. when a large number of random variables is involved). This methodology significantly reduces the computational cost with respect to the classical full polynomial chaos expansion (PCE) methodology. Notice however that when dealing with computationally-expensive deterministic models, the time cost remains important even with the use of the SPCE. In this paper, an efficient combined use of the SPCE methodology and the global sensitivity analysis (GSA) is proposed to solve such a problem. The proposed methodology is validated using a relatively non-expensive deterministic model.

INTRODUCTION

An efficient approach to deal with uncertainties propagation in case of high-dimensional problems (i.e. when a large number of random variables is involved) was recently presented by [START_REF] Blatman | An adaptive algorithm to build up sparse polynomial chaos expansions for stochastic finite element analysis[END_REF]. This approach is based on a Sparse Polynomial Chaos Expansion (SPCE) for the system response and leads to a reduced computational cost as compared to the classical Polynomial Chaos Expansion (PCE) methodology. Notice that both, the PCE and the SPCE methodologies, aim at replacing the original expensive deterministic model which may be an analytical model or a finite element/finite difference model by a meta-model. This allows one to calculate the system response (when performing a probabilistic analysis based on Monte Carlo Simulation MCS) using a simple analytical equation (e.g. [START_REF] Isukapalli | Stochastic response surface methods (SRSMs) for uncertainty propagation: Application to environmental and biological systems[END_REF][START_REF] Huang | Geotechnical probabilistic analysis by collocation-based stochastic response surface method: An Excel Add-in implementation[END_REF][START_REF] Mollon | Probabilistic analysis of pressurized tunnels against face stability using collocation-based stochastic response surface method[END_REF][START_REF] Mao | Probabilistic analysis and design of strip foundations resting on rocks obeying Hoek-Brown failure criterion[END_REF]. Notice however that when dealing with computationally-expensive deterministic models with a large number of random variables, the time cost remains important even with the use of the SPCE. Consequently, a method that can reduce once again the cost of the probabilistic analysis is needed. In this paper, an efficient combination between the SPCE methodology and the Global Sensitivity Analysis (GSA) is proposed to solve such a problem. In this method, a small SPCE order is firstly selected to approximate the system response by a meta-model. A GSA based on Sobol indices is then performed on this small SPCE order to determine the weight of each random variable in the variability of the system response. As a result, the variables with very small values of their Sobol indices (i.e. those that have a small weight in the variability of the system response) can be discarded. Consequently, a response which only depends on a smaller number of random variables is obtained. In other words, one obtains a response with an 'effective dimension'. This dimension is smaller than the initial dimension where the total number of random variables was considered. As it will be shown later, the use of a small SPCE order to perform the GSA is not a concern since higher SPCE orders lead to the same influential random variables. Once the 'effective dimension' was determined, a higher SPCE order that makes use of only the most influential random variables can be used. This significantly reduces the computation time. The use of a higher SPCE order is necessary in order to lead to an improved fit of the SPCE.

The proposed methodology is validated using a relatively non-expensive model which was extensively investigated by Al-Bittar andSoubra (2011, 2012). This model involves the computation of the ultimate bearing capacity of a strip footing resting on a weightless spatially varying (c, φ) soil where c is the soil cohesion and φ is the soil angle of internal friction. It should be noticed here that the random fields of c and φ are discretized into a finite number of random variables. This number is small for very large autocorrelation distances and significantly increases for small values of the autocorrelation distances.

The paper is organized as follows: The next two sections aim at briefly presenting both the sparse polynomial chaos expansion (SPCE) and the global sensitivity analysis (GSA). Then, the proposed efficient combination between the SPCE methodology and the GSA is presented. It is followed by the numerical results. The paper ends with a conclusion.

SPARSE POLYNOMIAL CHAOS EXPANSION (SPCE) METHODOLOGY

In this section, one first presents the polynomial chaos expansion (PCE) and then its extension, the sparse polynomial chaos expansion (SPCE). The Polynomial Chaos Expansion (PCE) methodology allows one to replace an expensive deterministic model which may be an analytical model or a finite element/finite difference numerical model by a meta-model. Thus, the system response may be calculated (when performing a probabilistic analysis based on Monte Carlo Simulation MCS) using a simple analytical equation. This equation is obtained by expanding the system response on a suitable basis which is a series of multivariate polynomials that are orthogonal with respect to the joint probability density function of the random variables. The PCE theory was originally formulated with standard Gaussian random variables and Hermite polynomials [START_REF] Spanos | Stochastic finite element expansion for random media[END_REF]. It was later extended to other types of random variables that use other types of polynomials [START_REF] Xiu | The Wiener-Askey polynomial chaos for stochastic differential equations[END_REF]. In this paper, standard normal random variables in conjunction with Hermite polynomials are used. The coefficients of the PCE may be efficiently computed using a non-intrusive technique where the deterministic calculations are done using for example an analytical model or a finite element/finite difference software treated as a black box. The most used non-intrusive method is the regression approach (e.g. [START_REF] Isukapalli | Stochastic response surface methods (SRSMs) for uncertainty propagation: Application to environmental and biological systems[END_REF][START_REF] Huang | Geotechnical probabilistic analysis by collocation-based stochastic response surface method: An Excel Add-in implementation[END_REF][START_REF] Blatman | An adaptive algorithm to build up sparse polynomial chaos expansions for stochastic finite element analysis[END_REF][START_REF] Mollon | Probabilistic analysis of pressurized tunnels against face stability using collocation-based stochastic response surface method[END_REF][START_REF] Mao | Probabilistic analysis and design of strip foundations resting on rocks obeying Hoek-Brown failure criterion[END_REF]. This method is used in the present work. The PCE methodology can be briefly described as follows:

For a deterministic model Γ with M random variables, the system response can be expressed by a PCE of order p fixed by the user as follows:

1 0 0 ( ) ( ) ( ) P PCE Y a a β β β β β β ξ ξ ξ ∞ - = = = Γ = Ψ ≅ Ψ   ( 1
)
where P is the number of terms retained in the truncation scheme,

{ } 1,...., i i M ξ ξ = =
is a vector of M independent standard random variables that represent the M random variables, a β are unknown coefficients to be computed and β Ψ are multivariate Hermite polynomials which are orthogonal with respect to the joint probability density function (PDF) of the standard normal random vector ξ. These multivariate Hermite polynomials can be obtained from the product of one-dimensional Hermite polynomials as follows:

1 ( ) = Ψ = ∏ i M i i H β α ξ (2)
where α i (i=1, …, M) are a sequence of M non-negative integers and (.)

i H α is the th i α one-
dimensional Hermite polynomial. The expressions of the one-dimensional Hermite polynomials are given in [START_REF] Spanos | Stochastic finite element expansion for random media[END_REF] among others.

In practice, the PCE with an infinite number of terms should be truncated by retaining only the multivariate polynomials β Ψ of degree less than or equal to p. For this purpose, the classical truncation scheme based on the determination of the first order norm 1 1 2011,2012). The SPCE methodology presented by [START_REF] Blatman | An adaptive algorithm to build up sparse polynomial chaos expansions for stochastic finite element analysis[END_REF] is an efficient alternative that can significantly reduce the number of calls of the deterministic model. In this methodology, [START_REF] Blatman | An adaptive algorithm to build up sparse polynomial chaos expansions for stochastic finite element analysis[END_REF] have shown that the number of significant terms in a PCE is relatively small since the multivariate polynomials β Ψ corresponding to high-order interaction (i.e. those resulting from the multiplication of the i H α with increasing α i values) are associated with very small values for the coefficients a β . Thus, a truncation strategy (called the hyperbolic truncation scheme) based on this observation was suggested by these authors. Within this strategy, the multivariate polynomials β Ψ corresponding to high-order interaction were penalized. This was performed by considering the hyperbolic truncation scheme which suggests that the q-norm should be less than or equal to the order p of the PCE. The q-norm is given by:

M i i α α = =  is used.
1 1 q M q i q i α α =   =      (3)
where q is a coefficient (0<q<1). In this formula, q can be chosen arbitrarily. [START_REF] Blatman | An adaptive algorithm to build up sparse polynomial chaos expansions for stochastic finite element analysis[END_REF] have shown that sufficient accuracy is obtained for 0.5 q ≥ . The proposed SPCE methodology leads to a sparse polynomial chaos expansion that contains a smaller number of unknown coefficients which can be calculated from a reduced number of calls of the deterministic model. This strategy was used in Al-Bittar andSoubra (2011, 2012) and will also be used in this paper to build up a SPCE of the system response. The iterative procedure suggested by [START_REF] Blatman | An adaptive algorithm to build up sparse polynomial chaos expansions for stochastic finite element analysis[END_REF] for building up a SPCE is detailed in [START_REF] Blatman | An adaptive algorithm to build up sparse polynomial chaos expansions for stochastic finite element analysis[END_REF] and Al-Bittar andSoubra (2011, 2012) and is not repeated herein. Once the coefficients a β have been computed, the statistical moments (mean, standard deviation, skewness, and kurtosis) can be calculated with no additional cost by generating a large number of realizations of the vector

{ } 1,...., i i M ξ ξ = =
and by computing the corresponding system responses using the obtained SPCE. The next subsection is devoted to the method used for the computation of the coefficients a β of the SPCE using the regression approach.

Computation of the SPCE coefficients by the regression approach

Consider a set of K realizations The computation of the SPCE coefficients using the regression approach is performed using the following equation:

( ) ( ) (
1 ( ) T T a η η η - = Γ 
(4) where the data matrix η is defined by:

( ) ( ), 1,..., , 0,..., 1 = Ψ = = - i i i K P β β η ξ β (5)
In order to ensure the numerical stability of the treated problem in Eq.( 4), the size K of the ED must be selected in such a way that the matrix 1 ( ) T η η -is well-conditioned. This implies that the rank of this matrix should be larger than or equal to the number of unknown coefficients. This test was systematically performed while solving the system of equations of the regression approach.

The quality of the output approximation via a SPCE closely depends on the SPCE order p. To ensure a good fit between the meta-model and the true deterministic model (i.e. to obtain the optimal SPCE order), the well-known coefficient of determination R 2 is generally used:

( ) ( ) ( ) ( ) ( ) ( ) 2 K i i SPCE i 1 2 2 K i i 1 1 K R 1 1 K 1 Γ ξ Γ ξ Γ ξ Γ = =   -   = -   -   -   Where (6) ( ) ( ) K i i 1 1 K Γ Γ ξ = =  (7)
The value 2 R 1 = indicates a perfect fit of the true model response Γ, whereas 2 R 0 =

indicates a nonlinear relationship between the true model Γ and the SPCE model SPCE

Γ

. The coefficient R 2 may be a biased estimate since it does not take into account the robustness of the meta-model (i.e. its capability of correctly predicting the model response at any point which does not belong to the experimental design). As a consequence, one makes use of a more reliable and rigorous coefficient of determination, namely the leave-one-out coefficient of determination [START_REF] Blatman | An adaptive algorithm to build up sparse polynomial chaos expansions for stochastic finite element analysis[END_REF]. This coefficient of determination consists in sequentially removing a point from the experiment design composed of K points. Let \i ξ Γ be the meta-model that has been built from the experiment design after removing the i th observation and let

i ( i ) ( i ) \i ( ) ( ) ξ Δ Γ ξ Γ ξ = -
be the predicted residual between the model evaluation at point

( i )
ξ and its prediction based on \i ξ Γ . The corresponding coefficient of determination is often

denoted by Q 2 : ( ) ( ) ( ) 2 K i i 1 2 2 K i i 1 1 K Q 1 1 K 1 Δ Γ ξ Γ = = = -   -   -   (8)
This coefficient will be used in the present paper to check the accuracy of the fit.

GLOBAL SENSITIVITY ANALYSIS (GSA)

Once the SPCE coefficients are determined, a global sensitivity analysis (GSA) based on Sobol indices can be easily performed. Notice that the first order Sobol index of a given random variable i ξ (i=1,…, M) gives the contribution of this variable in the variability of the system response. The first order Sobol index is given by [START_REF] Saltelli | Sensitivity Analysis[END_REF] and [START_REF] Sobol | Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates[END_REF] as follows:

( ) ( ) | ( ) i i V ar E Y S V ar Y ξ ξ     = (9)
where Y is the system response, ( )

| i E Y ξ
is the expectation of Y conditional on a fixed value of i ξ , and Var denotes the variance.

In the present paper, the system response Y is represented by a SPCE. Thus, by replacing Y in Eq.( 9) with the SPCE expression, one obtains the Sobol index formula as a function of the different terms of the SPCE [START_REF] Sudret | Global sensitivity analysis using polynomial chaos expansion[END_REF]. This formula is given by:

( ) ( ) 2 2 ( ) i I i PC a E S D β β β ξ ∈   Ψ     =  (10)
where a β are the obtained SPCE coefficients, β Ψ are the multivariate Hermite polynomials,

[ ]

. E is the expectation operator, and PC D is the variance of the response approximated by the SPCE. The response variance D PC is given by [START_REF] Sudret | Global sensitivity analysis using polynomial chaos expansion[END_REF] as follows:

( ) [ ] ( ) 1 2 2 2 , ,..., M PC I I I D a E β β β∈   = Ψ      (11)
Notice that the term ( )

2 E β   Ψ    
that appears in both Eq. ( 10) and Eq. ( 11) is given by Sudret (2008) as follows:

( )

2 1 ! M i i E β α = Ψ = ∏ (12)
where the α i are the same sequence of M non-negative integers { } 1 , ....., M α α used in Eq. ( 2).

Notice finally that I i in Eq. ( 10) denotes the set of indices β for which the corresponding β Ψ is only a function of the random variable i ξ (i.e. it only contains the variable i ξ ), and I i (i=1, …, M) regroup all the indices β for which the corresponding β Ψ is only a function of the random variable i ξ (i=1, …, M).

In order to illustrate the PCE theory and the global sensitivity analysis based on Sobol indices in a simple manner, an illustrative example of a PCE of order p=3 using only M=2 random variables ( 1 ξ and 2 ξ ) may be found in Al-Bittar and Soubra (2014).

EFFICIENT COMBINATION BETWEEN THE SPCE METHODOLOGY AND THE GLOBAL SENSITIVITY ANALYSIS

As mentioned previously, the time cost of the probabilistic analysis remains important even with the use of the SPCE when dealing with computationally-expensive deterministic models. Consequently, a procedure that can reduce once again this time cost is needed. An efficient combination between the SPCE methodology and the GSA is proposed in this section. The basic idea of this combination is that, for a given discretized random field, the obtained random variables do not have the same weight in the variability of the system response. The variables with a very small contribution in the variability of the system response can be discarded which significantly reduces the dimensionality of the treated problem. This allows one to perform a probabilistic analysis using a reduced Experiment Design (ED) and thus a smaller number of calls of the computationally-expensive deterministic model. The main challenge remains in detecting the most influential random variables in order to reduce the dimensionality of the problem. For this purpose, a procedure that makes use of both the SPCE and the GSA (denoted hereafter by SPCE/GSA) is proposed in this regard. The SPCE/GSA procedure can be summarized by the following steps:

Step 1: discretize the random field(s): This step was made in this paper using EOLE method and its extensions by Vořechovsky (2008). Let us consider N RF anisotropic non-Gaussian crosscorrelated random fields ( , )

NG i Z x y ( 1,..., RF i N =
) described by: (i) constant means and standard deviations (μ i , σ i ; 1,..., RF i N =

), (ii) non-Gaussian marginal cumulative distribution functions CDFs named G i ( 1,..., RF i N =

), (iii) a target cross-correlation matrix C NG and (iv) a common square exponential autocorrelation function NG Z ρ [(x, y), (x', y')] which gives the values of the correlation function between two arbitrary points (x, y) and (x', y'). This autocorrelation function is given as follows:

2 2 ' ' [( , ), ( ', ')] exp Z NG x y x x y y x y x y a a ρ       - -   = - -              ( 13
)
where a x and a y are the autocorrelation distances along x and y respectively. The Expansion Optimal Linear Estimation method (EOLE) and its extension by Vořechovsky ( 2008) to cover the case of correlated non-Gaussian random fields are used herein to generate the N RF random fields.

Notice that EOLE was first proposed by [START_REF] Li | Optimal discretization of random fields[END_REF] for the case of uncorrelated Gaussian fields, and then extended by Vořechovsky ( 2008) to cover the case of correlated non-Gaussian fields. In this method, one should first define a stochastic grid composed of q grid points (or nodes) { } 1 1 ( , ), ..., ( , ) q q

x y x y for which the values of the field are assembled in a vector { }

1 1 ( , ), ..., ( , )

q q Z x y Z x y χ =
. Then, one should determine the common correlation matrix for which each element ( )

, ; i j NG χ χ
Σ is calculated as follows:

( ) , ;

( , ), ( , )

Z i j NG NG i i j j x y x y χ χ ρ   Σ =   (14)
The common non-Gaussian autocorrelation matrix ; NG χ χ Σ and the target non-Gaussian crosscorrelation matrix C NG should be transformed into the Gaussian space using Nataf model [START_REF] Nataf | Détermination des distributions de probabilités dont les marges sont données[END_REF] since the discretization of the random fields using EOLE is done in the Gaussian space. As a result, one obtains N RF Gaussian autocorrelation matrices

; i χ χ Σ ( 1,..., RF i N =
), and a Gaussian cross-correlation matrix C that can be used to discretize the two random fields. The value i Z  of a random field obtained using this method is given by the following equation (cf. Al-Bittar andSoubra, 2011, 2012):

( ) ( , ); , 1 ( , ) . . 1, ..., Z x y j D N T i j i i i i i j R F i j Z x y µ i N χ κ σ φ λ = ≅ + Σ =   ( 15
)
where N RF is the number of random fields, N is the number of terms retained in the series expansion, ,

D i j

κ are N RF cross-correlated blocks of independent standard normal random variables obtained using the Gaussian cross-correlation matrix C between the N RF fields, ( , ; 1,..., Once the N RF Gaussian random fields are obtained, they should be transformed into the non-Gaussian space (in case of non-Gaussian random fields) by applying the following formula:

{ } 1 ( , ) ( , ) 1, ..., NG i i i R F Z x y G Z x y i N -   = Φ =     (16)
where (.) Φ is the standard normal cumulative density function (CDF). For more details about the EOLE method and its extensions to cover the case of cross-correlated non-Gaussian random fields, the reader may refer to Vořechovsky (2008) and Al-Bittar andSoubra (2011, 2012). After the discretization procedure, a random field is represented by N independent standard normal random variables. For the N RF random fields that have the same autocorrelation function, the total number of random variables is N T = N RF xN which can be relatively large especially for small values of the autocorrelation distances.

Step 2: use a preliminary small order of the sparse polynomial chaos expansion (e.g. p=2) to approximate the system response by a meta-model. The main reason for selecting a small SPCE order is the exploration of the most influential random variables (i.e. those that have a significant weight in the variability of the system response) using a small Experiment Design (ED). It should be emphasized here that the reduced number of the unknown SPCE coefficients related to the small value of the SPCE order leads to a significant decrease in the size of the experiment design, i.e. in the number of calls of the deterministic model.

Step 3: perform a GSA based on Sobol indices (using the obtained second order SPCE) to determine the weight of each random variable in the variability of the system response. The variables with very small values of their Sobol indices have no significant impact in the variability of the system response and can thus be discarded. Consequently, a response that only depends on a smaller number of random variables is obtained. In other words, one obtains a response with an 'effective dimension' N e that is smaller than the initial dimension where the total number N T of random variables was considered. It should be mentioned here that the small SPCE order (i.e. p=2) used firstly in step 2 to perform the GSA is sufficient to provide the weight of each random variable in the variability of the system response. This is because higher SPCE orders lead to the same influential random variables as will be seen later in the numerical results.

Step 4: use the same Experiment Design (ED) which was employed in step 2 but this time by only keeping the most influential random variables. By reducing the number of random variables from N T to N e (N e < N T ), one has the possibility to use a higher SPCE order (i.e. p>2). The use of a higher SPCE order is necessary to lead to an improved fit of the SPCE since the leave-one-out coefficient of determination Q 2 given in Eq. ( 8) increases when the SPCE order increases as it will be shown in the numerical results.

As a conclusion, the use of the SPCE/GSA procedure has the advantage of performing a good fit of the deterministic model with a reduced number of model evaluations as compared to the classical SPCE approach.

NUMERICAL RESULTS

The aim of this section is to validate the present SPCE/GSA procedure. For this purpose, a comparison between the results obtained by using the classical SPCE methodology and those given by the proposed SPCE/GSA procedure is presented. A computationally non expensive deterministic model was used for the validation. The problem used for the validation was presented in Al-Bittar andSoubra (2011, 2012). It aims at computing the ultimate bearing capacity of a strip footing resting on a c-φ spatially varying soil. The input parameters are similar to those considered in Al-Bittar andSoubra (2011, 2012). They are briefly presented in Table 1. For a more detailed description on these data, the reader may refer to Al-Bittar andSoubra (2011, 2012). It should be mentioned here that the friction angle φ and the dilation angle ψ were assumed to be perfectly correlated random fields since ψ =2φ/3. Thus, only the friction angle was modelled by a random field and the values of the dilation angle within the different elements of the deterministic mesh were deduced from the corresponding values of the soil angle of internal friction using the relationship ψ =2φ/3. The deterministic model is based on numerical simulations using FLAC 3D and it involves the case of a weightless soil. Thus, one obtains the soil bearing pressure due to only the soil cohesion; the contribution of the soil weight being neglected in the present paper. It should be mentioned here that when neglecting the soil weight γ, the computation time decreases from 10 to 5 min per simulation. This significantly reduces the computation time for the validation of the present SPCE/GSA procedure. As shown in Figure 1, the adopted soil domain considered in the analysis is 15m wide by 6m deep. For the boundary conditions, the horizontal movement on the vertical boundaries of the grid is restrained, while the base of the grid is not allowed to move in both the horizontal and the vertical directions.

Table 1. Deterministic and probabilistic input parameters

Lognormal

Reference case: The validation of the SPCE/GSA procedure is done for the illustrative case [a x =10m, a y =1m, r(c, φ) =-0.5] referred to hereafter as the reference case. For this configuration, the discretization of the two random fields c and φ has led to a total number of random variables N T equal to 24 (12 random variables for each random field as was shown in Al-Bittar and Soubra (2011, 2012)). By using the total number of random variables N T , Al-Bittar andSoubra (2011, 2012) have shown that a third order SPCE was sufficient to reach a target accuracy of 0.999. An ED involving 800 points was needed to solve the regression problem given in Eq. (4) (i.e. to obtain a wellconditioned regression problem for which the rank of the matrix 1 ( ) T η η -is larger than or equal to the number of unknown coefficients). By using the present SPCE/GSA procedure, a GSA was performed to detect the most influential random variables. Different SPCE orders (i.e. orders 2, 3, and 4) were considered in order to check if the SPCE order has an impact on the most influential random variables. Figure 2 depicts the values of Sobol indices for the 24 random variables, as given by SPCEs of orders 2, 3 and 4. The first 12 random variables [i.e. ξ i for i=1, …, 12] correspond to the cohesion random field and the last 12 random variables [i.e. ξ i for i=13, …, 24] are those corresponding to the friction angle random field. Figure 2 shows that whatever the SPCE order is, the two first random variables of both fields, (i.e. ξ 1 , ξ 2 , ξ 13 , ξ 14 ) are the most influential. For the two random fields, a very fast decay in the weight of the random variables is noticed with quasi negligible values beyond the first two random variables. In fact, the first two random variables of the two random fields, which correspond to the first two eigenmodes of both fields involve 95% of the response variability as may be seen from Table 2. This is logical since the system response (i.e. the ultimate bearing capacity) depends on the average distributions of the soil shear strength parameters c and φ over the soil domain. It is thus quite insensitive to the small-scale fluctuations of theses parameters around their average distributions.

( ) ( ) ( ) ( ) , 0 , 0 . 5 , 0 . 5 , 0 NG r c c r c C r c r ϕ ϕ ϕ ϕ   = = - =   = - =   Square exponential 2 2 ' ' exp NG Z x y x x y y a a ρ       - -   = - -             
Figure 2 clearly shows that the Sobol indices of the different random variables do not significantly change with the SPCE order. Thus, a second order SPCE is sufficient to identify the influential random variables (i.e. those that have a significant weight in the variability of the ultimate bearing capacity). The main advantage of a small SPCE order is that a small ED is sufficient to solve the regression problem. As shown in Table 3, 150 calls of the deterministic model are needed to solve the regression problem for a second order SPCE. This number attains 800 for a fourth order SPCE. This significant increase is because the number of unknown coefficients significantly increases from 29 to 144 when one chooses a fourth SPCE order instead of a second SPCE order. To choose the number of random variables which will be retained hereafter, the different random variables of the two random fields are firstly sorted in a descending order according to the values of their Sobol indices (cf. first and second columns in Table 4). A threshold of acceptance t a is then fixed as a percentage of the most influential (weighted) random variable. In the present paper, the most influential random variable is ξ 1 and it has a Sobol index equal to 0.5. Different values of the threshold were tested (cf. first line in Table 4). The random variables having a Sobol index smaller than the prescribed threshold t a are discarded. In this paper, a threshold of 2% of the Sobol index of the most weighed random variable is considered as sufficient; the corresponding retained random variables provide 98% of the total variance of the system response as may be seen from column 6 of Table 4. For this threshold, an 'effective dimension' N e =5 is obtained (i.e. 5 random variables are considered to be the most weighed). The 5 retained random variables will now be used with the already existing 150 model evaluations which were firstly employed to approximate the second order SPCE with the total number of random variables N T =24.

The reduction in the number of random variables from N T =24 to N e =5 provides the possibility to use higher SPCE orders (i.e. p>2) with the same ED (i.e. the 150 model evaluations). The use of a higher SPCE order is necessary to lead to an improved fit of the SPCE since the coefficient of determination Q 2 given in Eq. ( 8) increases when the SPCE order increases as shown in Table 5 for both the classical SPCE approach (using the total number of random variables N T =24) and the present SPCE/GSA procedure (where the effective dimension is equal to 5 (i.e. N e =5)). Using the SPCE/GSA procedure, an SPCE up to p=8 was reached using only 150 model evaluations. From Table 5, one can notice that with the use of the SPCE/GSA procedure, the Q 2 increases with the increase of the SPCE order and stabilizes beyond the order 5. This means that no improvement in the fit is obtained beyond this order. On the other hand, the value of Q 2 given by the present approach is smaller than the classical SPCE approach with a fourth order. This is because 19 random variables were discarded which slightly affect the goodness of the fit.

Figure 3 shows the PDF of the ultimate bearing capacity as obtained by both the classical SPCE approach (with the total number of random variables N T =24) and the proposed SPCE/GSA procedure (using only five random variables). Table 6 provides the corresponding statistical moments. This table also provides the coefficients of determination of the SPCEs. Notice that the results of the present SPCE/GSA approach are given in Table 6 for different values of the model evaluations (from 150 to 800). From this table, one can see that the coefficients of determination corresponding to the SPCE/GSA procedure are quasi constant with the increase in the number of model evaluations. This means that 150 model evaluations are sufficient and there is no need for more model evaluations to improve the accuracy of the fit. 3 and Table 6) that the first two statistical moments are well estimated with the present SPCE/GSA approach using the 150 model evaluations. However, the third and fourth statistical moments need more model evaluations (800 model evaluations) in order to converge to their reference values given by the SPCE approach (cf. Table 6). This demonstrates the efficiency of the present SPCE/GSA procedure to compute the first two statistical moments with a much reduced number of the model evaluations with respect to the classical SPCE approach.

As for the Sobol indices of the two random fields c and φ, Table 7 shows that the SPCE/GSA procedure with only 150 model evaluations gives the same results obtained by the classical SPCE approach using 800 model evaluations which demonstrates once again the efficiency of the present SPCE/GSA procedure.

CONCLUSIONS

An efficient combined use of the SPCE methodology and the global sensitivity analysis (GSA) has been proposed. The aim is to reduce the cost of the probabilistic analysis of computationally-expensive deterministic models. This methodology was validated in this paper using a relatively non-expensive deterministic model. The validation consists in comparing the results of both the classical SPCE methodology with the total number of random variables and the proposed combination between the SPCE and the GSA. Satisfactory results were obtained. A much smaller number of model evaluations was needed with the proposed methodology. The first two statistical moments and the Sobol indices have been well estimated with a small number of model evaluations. On the other hand, the third and fourth statistical moments need more model evaluations in order to converge to their reference values obtained using the classical SPCE. Since the present SPCE/GSA procedure was shown to be efficient for the probabilistic computation with a reduced calculation cost with respect to the classical SPCE approach, this approach may now be applied with confidence to costly deterministic models.

  vector between the value of the field i Z  at an arbitrary point (x, y) and its values at the different points of the stochastic mesh.

  Figure 1. The adopted soil domain

Figure 2 .

 2 Figure 2. Sobol indices for SPCEs of orders 2, 3 and 4 using the total number of eigenmodes ξ i (i=1, …, 24)

Figure 3 .

 3 Figure 3. PDF of the ultimate bearing capacity for both the classical SPCE with the total number of random variables N T =24 and the proposed SPCE/GSA procedure with only five random variables N e =5.

Table 2 .

 2 Sobol indices for the reference case where a x =10m, a y =1m, and r(c,φ)=-0.5

	Sobol index	ξ 1 0.50	ξ 2 0.17 0.002 ξ 3	ξ i (i=1, ..., 12) for the cohesion random field ξ 4 ξ 5 ξ 6 ξ 7 ξ 8 ξ 9 0.002 0.03 0.002 0.009 0.0002 0.0002 9 x10 -05 0.0002 7 x10 -05 ξ 10 ξ 11 ξ 12
	Sobol index	ξ 13 0.2	ξ i (i=13, ..., 24) for the friction angle random field ξ 16 ξ 17 ξ 18 ξ 19 ξ 20 ξ 21 0.08 0.001 0.0008 0.002 0.0005 0.0006 0.0003 0.0001 4 x10 -05 4 x10 -05 5 x10 -05 ξ 14 ξ 15 ξ 22 ξ 23 ξ 24

Table 3 .

 3 Number of unknown coefficients and model evaluations for different SPCE order

	SPCE order	2	3	4
	Number of unknown coefficients P	29	35	144
	Number of model evaluations	150	350	800

Table 4 .

 4 Sobol indices of the different random variables and the retained random variables for the different values of the threshold of acceptance

	Random variable	Sobol index	t a =0.5%×ξ 1 =0.0025	t a =1%×ξ 1 =0.005	t a =1.5%×ξ 1 =0.0075	t a =2%×ξ 1 =0.01	t a =2.5%×ξ 1 =0.0125	t a =3%×ξ 1 =0.015	t a =4%×ξ 1 =0.02	t a =5%×ξ 1 =0.025
	ξ 1 ξ 13 ξ 2 ξ 14 ξ 5 ξ 7 ξ 6 ξ 17 ξ 3 ξ 4 ξ 15 ξ 16 ξ 19 ξ 18 ξ 20 ξ 8 ξ 9 ξ 11 ξ 21 ξ 10 ξ 12 ξ 24 ξ 22 ξ 23	0.5 0.2 0.17 0.08 0.03 0.009 0.002 0.002 0.002 0.002 0.001 0.0008 0.0006 0.0005 0.0003 0.0002 0.0002 0.0002 0.0001 9.0 ×10 -5 7.0 ×10 -5 5.0 ×10 -5 4.0 ×10 -5 4.0 ×10 -5	0.5 0.2 0.17 0.08 0.03 0.009 0.002 0.002 0.002 0.002	0.5 0.2 0.17 0.08 0.03 0.009	0.5 0.2 0.17 0.08 0.03 0.009	0.5 0.2 0.17 0.08 0.03	0.5 0.2 0.17 0.08 0.03	0.5 0.2 0.17 0.08 0.03	0.5 0.2 0.17 0.08 0.03	0.5 0.2 0.17 0.08 0.03
	Sum of
	Sobol	1.001	0.997	0.989	0.989	0.98	0.98	0.98	0.98	0.98
	indices

Table 5 .

 5 SPCE using the total and the reduced number of random variables

		SPCE order	2	3	4	5	6	7	8
	Total number of	Coefficient of determination R 2	0.998 0.999 0.999	-	-	-	-
	random variables	Leave-one-out
	N T	coefficient of determination Q 2	0.824 0.932 0.994	-	-	-	-
	Reduced number	Coefficient of determination R 2	0.961 0.963 0.968 0.970 0.972 0.972 0.972
	of random	Leave-one-out
	variables N e	coefficient of determination Q 2	0.791 0.883 0.957 0.961 0.963 0.963 0.963
	It can be observed (see Figure

Table 6 .

 6 Coefficients of determination of the SPCEs and statistical moments of the ultimate bearing capacity as given by the classical SPCE approach and by the present SPCE/GSA procedure

								Number of model evaluations	Mean μ qult (kPa)	Standard deviation σ qult (kPa)	Skewness δ u (-)	Kurtosis κ u (-)	R 2	Q 2
	With the total	number of	random	variables	N T =24	800	658.2	93.57	0.287	0.163	0.999	0.995
								150	657.84	90.80	0.105	0.0129	0.972	0.957
	With the reduced	number of random	variables N e =5		200 250 300 400 500 600 700	658.98 659.90 659.73 660.05 659.50 659.75 659.50	91.53 92.10 92.15 90.95 90.81 90.99 90.85	0.168 0.188 0.202 0.291 0.296 0.272 0.280	0.0563 0.0630 0.0600 0.0500 0.0430 0.116 0.1637	0.972 0.964 0.962 0.969 0.970 0.968 0.968	0.951 0.956 0.963 0.960 0.963 0.963 0.963
								800	659.85	91.20	0.30	0.160	0.970	0.967

Table 7 .

 7 Sobol indices as computed from the classical SPCE approach and the present SPCE/GSA procedure.

			Number of model evaluations	i	S i (i=1, ..., 12)	i	S i (i=13, ..., 24)	( ) S c	12 1 i = = 	i S	S	( ) ϕ	24 13 = =  i	i S
	With the total number of random	T =24 variables N	800	1 2 3 4 5 6 7 8 9 10 11 12	0.5 0.17 0.002 0.002 0.03 0.002 0.009 0.0002 0.0002 9.0 ×10 -5 0.0002 7.0 ×10 -5	13 14 15 16 17 18 19 20 21 22 23 24	0.2 0.08 0.001 0.0008 0.002 0.0005 0.0006 0.0003 0.0001 4.0 ×10 -5 4.0 ×10 -5 5.0 ×10 -5	0.715			0.285
			Number of
			model evaluations	i S i (i=1, 2, 3) i	S i (i=4, 5)


Françoise Masson
email: f.masson@crossdata.it

Elias Cueto
email: ecueto@unizar.es

David Gonzalez
email: gonzal@unizar.es

Francisco Chinesta
email: francisco.chinesta@ec-nantes.fr

DYNAMIC DATA-DRIVEN PARAMETRIC DYNAMICAL SYSTEMS

Keywords: Dynamic Data-Driven Application Systems, Model reduction, PGD

Dynamic Data-Driven Application Systems -DDDAS -appear as a new paradigm in the field of applied sciences and engineering, and in particular in simulation-based engineering sciences. By DDDAS we mean a set of techniques that allow the linkage of simulation tools with measurement devices for real-time control of systems and processes. DDDAS entails the ability to dynamically incorporate additional data into an executing application, and in reverse, the ability of an application to dynamically steer the measurement process. DDDAS needs for accurate and fast simulation tools making use if possible of off-line computations for limiting as much as possible the on-line computations. We could define efficient solvers by introducing all the sources of variability as extra-coordinates in order to solve the model off-line only once to obtain its most general solution to be then considered in on-line purposes. In this work we are evaluating one such strategy in both linear and non-linear dynamical systems as the ones usually encountered in control, epidemiology models, physiological systems (diagnostic of diabetes …). Such models result defined in highly multidimensional spaces suffering the so-called curse of dimensionality that could be efficiently circumvented by applying a Proper Generalized Decomposition that constitutes a work in progress.

INTRODUCTION

Traditionally, Simulation-based Engineering Sciences -SBES-relied on the use of static data inputs to perform the simulations. These data could be parameters of the model(s) or boundary conditions, outputs at different time instants, etc., traditionally obtained through experiments. The word static is intended to mean here that these data could not be modified during the simulation.

A new paradigm in the field of Applied Sciences and Engineering has emerged in the last decade. Dynamic Data-Driven Application Systems (DDDAS) constitute nowadays one of the most challenging applications of simulation-based Engineering Sciences. By DDDAS we mean a set of techniques that allow the linkage of simulation tools with measurement devices for real-time control of simulations. As defined by the U.S. National Science Foundation, "DDDAS entails the ability to dynamically incorporate additional data into an executing application, and in reverse, the ability of an application to dynamically steer the measurement process" [START_REF]DDDAS Workshop[END_REF] The term Dynamic Data-Driven Application System was coined by F. Dalema in a NSF workshop on the topic in 2000. The document that initially put forth this initiative stated that DDDAS constitute "application simulations that can dynamically accept and respond to 'online' field data and measurements and/or control such measurements. This synergistic and symbiotic feedback control loop among applications, simulations, and measurements is a novel technical direction that can open new domains in the capabilities of simulations with a high potential payoff, and create applications with new and enhanced capabilities. It has the potential to transform the way science and engineering are done, and induces a major beneficial impact in the way many functions in our society are conducted, such as manufacturing, commerce, transportation, hazard prediction/management, and medicine, to name a few" [START_REF] Darema | Engineering/Scientific and Commercial applications: differences, similarities, and future evolution[END_REF] A DDDAS includes different constituent blocks:

1. A set of (possibly) heterogeneous simulation models. 2. A system to handle data obtained from both static and dynamic sources. 3. Algorithms to efficiently predict system behavior by solving the models under the restrictions set by the data. 4. Software infrastructure to integrate the data, model predictions, control algorithms, etc.

Almost a decade after the establishment of the concept, the importance of the challenge is now better appreciated. As can be noticed, it deals with very different and transversal disciplines: from simulation techniques, numerical issues, control, modeling, software engineering, data management and telecommunications, among others. In a nutshell, the characteristics of a DDDAS could be depicted as in Fig. 1. Figure 1 identifies three different blocks of interactions:

(i) the one between human systems and the simulation, (ii) the simulation interaction with the physical system and (iii) the simulation and the hardware/ data infrastructure. Physical systems operate at very different time scales. In the picture (adapted from the NSF workshop on DDDAS 2000 final report) we have included some representative values: from 10 -20

Hz for cosmological systems to 10 20 Hz for problems at the atomic scales. Humans, however, can be considered as a system operating at rates from 3 Hz to 500Hz in haptic 1 devices for instance to transmit realistic touch sensations. A crucial aspect of DDDAS is that of real-time simulation. This means that the simulations must run at the same time (or faster) than the data are collected. While this is not always true (as in weather forecasting, for instance, where collected data are not usually incorporated to the simulations), most applications require different forms of real-time simulations. In haptic surgery simulators, for instance, the simulation result, i.e., forces acting on the surgical tool, must be translated to the peripheral device at a rate of 500 Hz, which is the frequency of the free hand oscillation. In other applications, such as some manufacturing processes, the time scales are much bigger, and therefore real-time simulations can last for seconds or minutes.

As can be noticed from the introduction above, DDDAS can revolutionize the way in which simulation will be done in the next decades. No longer a single run of a simulation will be considered as a means of validate a design on the basis of a static data set.

The importance of DDDAS in the forthcoming decades can be noticed from the NSF Blue Ribbon Panel on Simulation Based Engineering Sciences report, that in 2006 included DDDAS as one of the five core issues or challenges in the field for the next decade (together with multiscale simulation, model validation and verification, handling large data and visualization). This panel concluded that "Dynamic data-driven application systems will rewrite the book on the validation and verification of computer predictions" and that "research is needed to effectively use and integrate data-intensive computing systems, ubiquitous sensors and high-resolution detectors, imaging devices, and other data-gathering storage and distribution devices, and to develop methodologies and theoretical frameworks for their integration into simulation systems" [START_REF] Oden | Simulation-based Engineering Science: Revolutionizing Engineering Science through simulation[END_REF]. Moreover, the NSF believes that "… The DDDAS community needs to reach a critical mass both in terms of numbers of investigators, and in terms of the depth, breadth and maturity of constituent technologies..." [START_REF]DDDAS Workshop 2006[END_REF]. While research on DDDAS should involve applications, mathematical and statistical algorithms, measurement systems, and computer systems software methods, our work focuses on the development of mathematical and statistical algorithms for the simulation within the framework of such a system. In brief, we intend to incorporate a new generation of simulation techniques into the field, allowing to perform faster simulations, able to cope with uncertainty, multiscale phenomena, inverse problems and many other features that will be discussed. This new generation of simulation techniques has received the name of Proper Generalized Decompositions -PGD-and has received an increasing level of attention by the SBES community [START_REF] Chinesta | Recent advances in the use of the Proper Generalized Decomposition for solving multidimensional models[END_REF]. PGD was initially introduced for addressing multidimensional models encountered in science and engineering and was then extended to address general computational mechanics models.

DYNAMIC DATA-DRIVEN PARAMETRIC DYNAMICAL SYSTEMS

In this work we propose a dynamic parameter estimation procedure making use of a general solution obtained off-line. The first results allow us to conclude on the potentiality of such procedure to be employed in DDDAS.

The general off-line solution is obtained by solving a linear or non-linear parametric dynamical system by introducing the initial condition and all the evolving or uncertain model parameters as extra-coordinates. The resulting multidimensional model is then solved only once, and from this solution we can evaluate the field of interest for any value of the model parameters (that can evolve in time) and the initial condition. In this framework, the information coming from measurements can be incorporated in real time to the on-line simulation.

The introduction of extra-coordinates was the price to be paid in order to compute a solution general enough to be used on-line without the necessity of further integrations of the dynamical system. On-line tasks only involve the post-processing of such general solution. Models involving many fields and many parameters limit the use of standard discretization techniques because the resulting model multidimensionality. The necessity of using fine enough meshes for representing the solution makes impossible its integration and even the data storage. A possibility for circumventing the difficulties relates to the curse of dimensionality lies in the use of the proper generalized decomposition proposed in some of our former works [START_REF] Chinesta | Recent advances in the use of the Proper Generalized Decomposition for solving multidimensional models[END_REF].

Figure 1 :

 1 Figure 1: Structure of a Dynamic Data-Driven Application System and typical rate of communication between blocks

Haptic devices are those able to translate to the user force reactions, such as, for instance, those used in surgical training.