Real-time data assimilation and control on mechanical systems under uncertainties

This research work deals with the implementation of so-called Dynamic Data-Driven Application Systems (DDDAS) in structural mechanics activities. It aims at designing a real-time numerical feedback loop between a physical system of interest and its numerical simulator, so that (i) the simulation model is dynamically updated from sequential and in situ observations on the system; (ii) the system is appropriately driven and controlled in service using predictions given by the simulator. In order to build such a feedback loop and take various uncertainties into account, a suitable stochastic framework is considered for both data assimilation and control, with the propagation of these uncertainties from model updating up to command synthesis by using a specific and attractive sampling technique. Furthermore, reduced order modeling based on the Proper Generalized Decomposition (PGD) technique is used all along the process in order to reach the real-time constraint. This permits fast multi-query evaluations and predictions, by means of the parametrized physics-based model, in the online phase of the feedback loop. The control of a fusion welding process under various scenarios is considered to illustrate the proposed methodology and to assess the performance of the associated numerical architecture.


Introduction
The continuous interaction between physical systems and high-fidelity simulation tools (i.e. virtual twins) has become a key enabler for industry as well as an appealing research topic along the last decade (see for instance [11]). This is at the heart of the Dynamic Data Driven Application System (DDDAS) concept [12], in which a simulation model is used to make decisions and drive an evolving physical system, and is in the same time fed by data collected on this system in order to update parameters and ensure the continual consistency between numerical predictions and physical reality. In other words, the DDDAS concept aims at building a numerical feedback loop between the physical system and its simulator, with on-the-fly data assimilation and control (Fig. 1). Nevertheless, there are two main numerical challenges in the implementation of such a loop for structural mechanics applications. On the one hand, the dialog between numerical models and physical systems is in practice subject to several sources of uncertainty, including measurement noise, modeling errors, or variabilities in the system properties and environment. On the other hand, a relevant feedback loop requires effective numerical methods such that real-time computations and interactions can be performed.
The paper presents a general strategy, addressing the two previous challenges, for the design of an effective numerical feedback loop between a physical system and its simulator. It considers a stochastic framework for sequential data assimilation and control, that uses Bayesian inference for model updating from in situ data as well as uncertainty propagation to make predictions from the model and synthesize control laws. Such a framework considers parameters to be inferred as random variables, and it naturally takes all uncertainty sources into account [2,6,17,22,30,31]. The proposed strategy also leans on two ingredients which permit to achieve the real-time constraint. First, Transport Map sampling [13] is used as an alternative to Markov Chain Monte-Carlo (MCMC) [14,25] or Sequential Monte-Carlo [1] techniques in order to perform fast Bayesian inference with convenient sampling of multi-dimensional posterior densities and associated adaptive strategies. The Transport Map technique consists in building a deterministic polynomial mapping between the posterior probability measure of interest and a simple reference measure (e.g. Gaussian distribution) [21,23,29]. It thus permits an automatic exploration, from the constructed mapping, of the multi-parametric stochastic space in order to effectively derive useful information such as means, standard deviations, maxima, or marginals on model parameters. Such pieces of information can then be propagated to model outputs in order to quantify uncertainty, synthesize the appropriate command in a stochastic context, and thus make safe decision on the evolving system. Second, model reduction by means of the Proper Generalized Decomposition (PGD) technique [9,10] is introduced in order to reduce the computational effort for the evaluation of multi-parametric numerical models, and therefore further speed up the overall process. The PGD approximation builds a modal representation of the multi-parametric model solution with separated variables and explicit dependency on model parameters. This representation is computed in an offline phase with controlled accuracy [8] before being evaluated at low cost in the online phase. It is shown in the paper that the PGD technique (i) facilitates the computation of the likelihood function involved in the Bayesian inference framework [3,26]; (ii) can be effectively coupled with Transport Map sampling for the calculation of the maps, as it directly provides information on solution derivatives [27,28]; (iii) is a particularly effective tool for performing uncertainty propagation through the forward model as well as command law synthesis. A particular focus is made here on the latter point dealing with effective command in a stochastic framework; this has been investigated in very few works of the literature, even though it is a major aspect of the DDDAS procedure. The dynamic command synthesis we propose, using advantages of Transport Map sampling and PGD model reduction, is the main novelty of the paper. It permits the construction and implementation of the full DDDAS feedback loop.
The constructed feedback loop is here illustrated in the context of a fusion welding process. It involves a simplified welding model introduced in [16] (and described in Fig. 2), which is supposed to be an accurate enough representation of the physical phenomena of interest.
In this two-dimensional model, two metal plates are welded by a heat source whose center is moving along the geometry. The problem unknown is the dimensionless temperature field T in the space domain and over the time domain I; T = 0 when the temperature is equal to the room temperature, and T = 1 when the temperature is equal to the melting temperature of the material. On the right-hand side boundary D (see Fig. 2), the temperature is supposed to be equal to the room temperature (T = 0). The other boundaries are supposed to be insulated. To solve the problem, the system of coordinates is made moving at the same speed as the heat source. Thus, the model problem is described by the following heat equation with convective term: where v = [Pe; 0] is the advection velocity, Pe = v · L c /κ is the Peclet number (L c being the characteristic length of the problem), and κ is the thermal diffusivity of the material. The volume heat source term s is defined by the following Gaussian repartition in the space domain: where coordinates (x c , y c ) represent the location of the heat source center, u is the magnitude, and σ is a scalar parameter that drives the source expansion. From the integration of (1) over , the weak formulation in space of the problem is of the form: find T ∈ T such that with: The functional space T is the Bochner space L 2 (I; S) S ⊗I, with S = H 1 0| D the Sobolev space of H 1 functions on satisfying homogeneous Dirichlet boundary conditions on D , and I = L 2 (I) the Lebesgue space.
The model parameters to be updated from indirect noisy data are p = {σ , Pe}, which are respectively related to the spatial spreading and speed of the heat source as illustrated in Fig. 3. They may be varying over the time domain. Data consist in the measurement of temperatures T 1 and T 2 at two points in (see Fig. 2). From these data assimilated sequentially in time, the purpose is twofold: (i) to dynamically update the model parameters p; (ii) to control from the updated model the temperature T 3 at another point in , which is the output of interest assumed to be unreachable by direct measurement, and perform corrections on the welding process if necessary. The control variable is the magnitude u of the heat source, that is supposed to be piecewise constant in time as illustrated in Fig. 3.
The paper outline is as follows: in "Reduced order modeling using PGD" section, the PGD model reduction applied to the above reference model is detailed. It is then employed in association with Bayesian inference and Transport Map sampling for fast data assimilation and model updating in "Real-time data assimilation with Bayesian inference and Transport Map sampling" section. All these tools are beneficially reused for on-the-fly command synthesis and system control in "Real-time control" section. Several numerical experiments are reported in "Results and discussion" section, which show the interest and performance of the proposed feedback loop by considering various welding scenarios. Sequential data assimilation, uncertainty propagation up to the output of interest, and real-time control of the welding process are illustrated for each of these scenarios. Eventually, conclusions and prospects are drawn in "Conclusions" section.

Reduced order modeling using PGD
Due to the increasing number of high-dimensional approximation problems, which naturally arise in many situations such as optimization or uncertainty quantification, model reduction techniques have been the object of a growing interest and are now a mature technology [19,24]. Tensor methods are among the most prominent tools for the construction of model reduction techniques as in many practical applications, the approximation of high-dimensional solutions of Partial Differential Equations (PDEs) is made computationally tractable by using low-rank tensor formats. In particular, an appealing technique based on a canonical format and referred to as Proper Generalized Decomposition (PGD) was introduced and successfully used in many applications of computational mechanics dealing with multiparametric problems [5,7,9,10,15,18,20]. Contrary to POD, the PGD approximation does not require any knowledge on the solution, and it operates in an iterative strategy in which basis functions (or modes) are computed from scratch by solving eigenvalue problems. In the classical PGD framework, the reduced model is built directly from the weak formulation (here (3)) of the considered PDE, integrated over the parametric space. The approximate reduced solution T m at order m is then is then searched in a in a separated form with respect to space, time, and model parameters p = {p 1 , p 2 , . . . , p d } seen as extra-coordinates [10]: The computation of the PGD modal representation is performed in an offline phase by using an iterative method [10], before being evaluated in an online phase at any space-time location and any parameter value from products and sums of one-parameter functions. For the multi-parametric problem of interest, the construction of the PGD solution is detailed in [26]. It reads: Considering a heat source term with u = 1, the first four PGD modes are represented in

Basics on Bayesian inference
The purpose of Bayesian inference is to characterize the posterior probability density function (pdf) π(p|d obs ) of some model parameters p given some indirect and noisy observations d obs . In this context, the Bayesian formulation of the inverse problem reads [17]: π(p|d obs ) = 1 C π(d obs |p).π 0 (p) where π 0 (p) is the prior pdf, related to the a priori knowledge on the parameters before the consideration of data d obs , π(d obs |p) is the likelihood function that corresponds to the probability for the model M to predict observations d obs given values of the parameters p, and C = π(d obs |p) · π(p)dp is a normalization constant. No assumption is made on the probability densities (prior, measurement noise) or on the linearity of the model. We consider here the classical case of an additive measurement noise with density π meas . We also consider that there is no modeling error, even though such an error source could be easily taken into account in the Bayesian inference framework (provided quantitative information on this error source is available). The likelihood function thus reads: Furthermore, when considering sequential assimilation of measurements d obs Once the PGD approximation T m (x, t, p) is built (see "Reduced order modeling using PGD" section), an explicit formulation of the non-normalized posterior density can be derived. Indeed, owing to the observation operator O, the output d m (p, t) = O (T m (x, t, p)) can be easily computed for any value of the parameter set p. The nonnormalized posterior density π thus reads: From the expression of π(p|d obs ) (or π(p|d obs 1 , . . . , d obs i )), stochastic features such as means, variances, or first-order marginals on parameters may be computed. These quantities are based on large dimension integrals, and classical Monte-Carlo integration-based techniques such as Markov Chain Monte-Carlo (MCMC) require in practice to sample the posterior density a large number of times. This multiquery procedure is much time consuming and incompatible with fast computations; we thus deal with an alternative approach in the following section.

Transport Map sampling
The principle of the Transport Map strategy is to build a deterministic mapping M between a reference probability measure ν ρ and a target measure ν π . The purpose is to find the change of variables such that:

Fig. 7 Illustration of the Transport Map principle for sampling a target density
In this framework, samples drawn according to the reference density are transported to become samples drawn according to the target density ( Fig. 7). For the considered inference methodology, the target density corresponds to the posterior density π(p|d obs ) derived from the Bayesian formulation, while a standard normal Gaussian density may be chosen as the reference density; for more details, we refer to [29] with effective computation tools (see http://transportmaps.mit.edu). From the reference density ρ, the purpose is thus to build the map M : R d → R d such that: where denotes the push forward operator. Once the map M is found, it can be used for sampling purposes by transporting samples drawn from ρ to samples drawn from π. Similarly, Gaussian quadrature (ω i , for π. Such a (deterministic) numerical integration with quadrature rule from the reference Gaussian density is therefore a technique of choice used in the present work for the calculation of statistics, marginals, or any other information from the posterior pdf. Maps M are searched among Knothe-Rosenblatt rearrangements (i.e lower triangular and monotonic maps). This particular choice of structure is motivated by the following properties (see [4,21,29] for all details): • Uniqueness and existence under mild conditions on ν π and ν ρ ; • Easily invertible map and Jacobian ∇M simple to evaluate; • Optimality regarding the weighted quadratic cost; • Monotonicity essentially one-dimensional (∂ p k M k > 0).
The maps M are therefore parametrized as: Functions c and e are chosen as Hermite polynomials with coefficients a c et a e . This integrated squared parametrization is a classical choice that automatically ensures the monotonicity of the map, and using Hermite polynomials leads to an integration that can be performed analytically. With this parametrization, the optimal map M is found by minimizing the following Kullback-Leibler (K-L) divergence: that quantifies the difference between the two distributions ν π and M ν ρ . Still using a Gaussian quadrature rule (ω i , p i ) N i=1 over the reference probability space associated with ρ, the minimization problem reads: where π is the non-normalized version of the target density. This minimization problem is fully deterministic and may be solved using classical algorithms (such as BFGS) using gradient or Hessian information on the density π (p). It is important to notice that the reduced PGD representation (6) of the solution is highly beneficial to solve (15). Partial derivatives of the model with respect to parameters p can indeed be easily computed as: and stored in the offline phase. Thanks to the separated representation of the PGD, crossderivatives are computed by combination of univariate modes derivatives. As a result, the use of PGD also speeds up the computation of transport maps. The quality of the approximation M ν ρ of the measure ν π can be estimated by the convergence criterion σ (variance diagnostic) defined in [29] as: The numerical cost for computing this criterion is very low as the integration is performed using the reference density and with the same quadrature rule as the one used in the computation of the K-L divergence. Therefore, an adaptive strategy regarding the order of the map can be used to derive an automatic algorithm that guarantees the quality of the approximation M ν ρ .
In the case of sequential inference, the Transport Map method exploits the Markov structure of the posterior density (9). Indeed, instead of being fully computed, the map between the reference density ρ and the posterior density at time t i is obtained by composition of low-order maps (see Fig. 8): Therefore, at each assimilation step t i , only the last map component M i is computed between ρ and the density π * i defined as: which leads to a process with almost constant CPU effort.

Real-time control
In addition to the mean, maximum a posteriori (MAP), or other estimates on model parameters, another major post-processing in the DDDAS feedback loop is the prediction of some quantities of interest from the model, such as the temperature T 3 at remote point x 3 in the present context (see Fig. 2). Once parameters p (σ and Pe here) are inferred in a probabilistic way at each assimilation time point t i (1 ≤ i ≤ N t ), it is indeed valuable to propagate uncertainties a posteriori in order to know their impact on the output of interest T 3 during the process, and consequently to assess the welding quality. As the PGD model gives an explicit prediction of the temperature field over the whole space-time-parametric domain, the output T 3 can be easily computed for all values of the parameter samples and at each physical time point τ j , j ∈ {1, . . . , N τ }. For a given physical time point τ j , the pdf π(T 3|τ j |p, t i ) of the value of the temperature T 3 knowing uncertainties on the parameter set p from data assimilation up to time point t i can thus be computed in real-time and used to determine if the plates are correctly welded and with which confidence. In practice, this computation may be performed for all physical time points τ j ≥ t i , and the density π(T 3|τ j |p, t i ) is characterized by a (Gaussian) quadrature rule using the Transport Map method. With this knowledge, a stochastic computation of the predicted temperature evolution can be obtained, and the control of the welding process from the numerical model can be performed. We detail below the procedure to dynamically determine the value of the control variable u (magnitude of the heat source) in the case where the welding objective is to satisfy a sufficient welding depth. The quantity of interest is then the maximal value of the temperature T 3 obtained at final time τ * , which is an indicator of the welding quality. When T 3|τ * ≥ 1, the welding depth is supposed to be sufficient. Other welding objectives will be considered in "Results and discussion" section, associated with similar strategies for command synthesis. Due to the stochastic framework which is employed, the quantity of interest is actually a random variable with pdf π(T 3|τ * |p, t i ) evolving at each data assimilation time t i . The proposed quantity q to monitor is: where Q is an operator defined in the stochastic space. This way, setting the objective q obj = 1 ensures that the temperature T 3|τ * is larger than the melting temperature with a confidence of 99%, and using the minimal energy (no overheating).
Using the PGD solution computed in "Reduced order modeling using PGD" section for a unit magnitude of the heat source (u = 1) and zero initial conditions, the predicted (stochastic) maximal value T 3 for a given constant magnitude u and for fixed pdfs of p reads: so that q = u · Q (T m (x 3 , τ * , p)) can be obtained in a straightforward manner. This way, setting the source magnitude u to u 0 = q obj /Q (T m (x 3 , τ * , p)) would enable to reach the welding objective. Nevertheless, in practice the pdfs on parameters p are updated at each assimilation time point t i , based on additional experimental information, so that the value of u needs to be tuned with time accordingly. In order to do so, the control variable u(t) is made piecewise constant in time, under the form: where H is the Heaviside function, u 0 is the initial command on the source magnitude (defined from the prior pdfs on p), and δu i is the correction to the current command at each assimilation time t i . Using the linearity of the problem with respect to the loading, a PGD solution associated with the command is made of a series of PGD solutions translated in time; it reads: Therefore, after each assimilation time point t i , the new prediction of the quantity of interest T 3|τ * can be easily obtained from PGD: where T pred,[0,i−1] 3|τ * (p) = u 0 · T m (x 3 , τ * , p) + i−1 n=1 δu n · T m (x 3 , τ * − t n , p) is the prediction on T 3|τ * considering the history of the control variable u(t) until time t i . Consequently, the correction δu i is defined such that Q(T 3|τ * ) = q obj , using (24) and considering the current pdfs of the parameter set p (i.e. those obtained after the last Bayesian data assimilation at time t i ).

Results and discussion
We now implement the DDDAS procedure proposed in "Methods" section on the model problem defined in "Introduction" section. We investigate three test cases involving different welding scenarios, in order to illustrate the flexibility of the approach and show its performance. For all scenarios, two temperature data T obs 1 and T obs 2 are assimilated at each assimilation time point t i in order to refine the knowledge on parameters σ and Pe, and further predict the value of the quantity of interest for control purpose. Without any limitation, we assume that assimilation time points t i , i ∈ {1, . . . , N t }, coincide with discretization time points τ j .

Case 1: control of the welding depth with constant physical process parameters
In this first test case, the control objective is the one mentioned in "Real-time control" section, that is Q(T 3|τ * ) = 1, with Q the operator defined in (20) and τ * = 45. This ensures that the temperature T 3 at final time τ * is larger than the melting temperature with a confidence of 99%, while using the minimal source energy. We use synthetic data, measurements being simulated using the PGD model with reference parameter values (σ ref = 0.4, Pe ref = −60) that are supposed to be constant in time in this section. An independent random normal noise is added with zero mean and standard deviations σ meas 1 = 0.01925 and σ meas 2 = 0.01245. Figure 9 shows the model outputs T 1 and T 2 at each time step as well as the perturbed outputs which provide the measurements used for the considered example, in the case where the control on the system is not activated (i.e. u = 1). When this control is implemented (see "On-the-fly control of the welding process" section), synthetic data are generated by taking into account the applied control law.
The goal of the test case is to perform a detailed analysis of the proposed DDDAS approach, in terms of dynamical model updating, uncertainty propagation on the quantity of interest, and on-the-fly command synthesis.

Dynamical updating of model parameters
The prior density on the parameters (σ , Pe) is chosen as the product of two independent Gaussian densities with means (μ σ = 0.4, μ Pe = −60) and variances (σ 2 σ = 0.003, σ 2 Pe = 7). The Transport Map strategy detailed in "Real-time data assimilation with Bayesian inference and Transport Map sampling" section and coupled with PGD is then applied for sequential data assimilation, assuming for the moment a constant magnitude u = 1 of the heat source. The solution of the heat equation (1) is used in its PGD form and derivatives of the approximate solution T m with respect to the parameters to be inferred are computed in order to derive the transport maps (i.e. successive maps M 1 , . . . , M N t ) effectively. In Table 1 we represent the computation time required to compute the transport maps at each assimilation step. We compare computation times when different information on derivative orders is provided to the minimization algorithm. With order 0, the minimization problem (15) is solved using a BFGS algorithm where the gradient is computed numerically. With order 1, the minimization is also performed using a BFGS algorithm but with the gradient given explicitly with respect to the PGD modes derivatives. With order 2, a conjugate gradient algorithm is used with an explicit formulation of both gradient and Hessian. The stopping criterion is a tolerance of 10 −3 on the variance diagnostic (17), and the complexity of the maps (order of the Hermite polynomials) is increased until this tolerance is fulfilled. It appears that the first assimilation step is the most expensive as the complexity of the transformation between the reference and the first posterior density is large (a 4th order map is required to fulfill the variance diagnostic criterion). The other transformations computed at other assimilation time steps are much less expensive (time less than 1 s) as they are built between intermediate posteriors which slightly differ at each step and can thus be easily represented by a linear (i.e. first order) transformation. The speed-up for the first iteration is about 5.5 between zeroth-order information and first-order information.
Between the first-order information and the second-order information, the speed-up is about 1.34. For the other time steps, the speed-up is very small as the computed map is very simple. We observe that using gradient and Hessian information to solve the minimization problem related to the computation of the transport maps leads to low computation times.
In Fig. 10, information on the computation cost over the time steps and using both gradient and Hessian information (order 2 information) is provided: Fig. 10a shows the computation time to build each map M i , i ∈ {1, . . . , N t }, while the cost in terms of model evaluations to compute each map is displayed in Fig. 10b. A level 10 Gauss-Hermite quadrature is used. From the second step to the final step, we observe that the computation time slowly increases (Fig. 10a) while the evaluation cost slowly decreases (Fig. 10b). This is due to the fact that the evaluation of the composition of maps grows with the number of steps. One way to circumvent this issue would consist in performing regression on the map composition. Figures 11 and 12 represent the marginals at each time step and for both parameters σ and Pe, respectively. The color map informs on the probability density function values.  During the iterations over the time steps, we observe that marginals become thinner with larger maximal pdf values giving more confidence on the parameters estimation. We also observe that the parameter σ is less sensitive than the parameter Pe regarding the inference process.

Uncertainty propagation on the quantity of interest
Still assuming a constant magnitude u = 1 of the heat source, uncertainty propagation is performed in real-time in order to predict the evolution of the temperature T 3 (in terms of pdf) in the region of interest. Knowing the uncertainties on the parameters, the goal is to predict at each assimilation time point the evolution of the temperature T 3 during the next physical time steps. This is easily done owing to the PGD model, as the temperature field is then globally and explicitly known over the time domain and with respect to the values of σ and Pe. The computation is performed after each assimilation time point t i and for all the physical time points τ j ≥ t i .  Other graphs (Fig. 13b-d) show the refinement of the prediction with the improvement on the parameters uncertainty knowledge. The current measurement assimilation step is indicated by the vertical cursor. On the right of the cursor τ = t i , the graphs represent the prediction of the temperature T 3 from the model after the assimilation of the measurements T obs,1:i 1 and T obs,1:i 2 . On the left of the cursor, each slice [t j−1 , t j ] (j ≤ i) represents the prediction made at the assimilation time t j (the predictions of the temperature T 3 for physical time steps anterior to the assimilation time step t i are not updated). Figure 14 shows the convergence of the prediction on the quantity of interest T 3|τ * at the steady state regime (τ * = 45) with respect to the assimilation steps. We observe that, as foreseen, more confidence is given to this output along the real-time data assimilation process.

On-the-fly control of the welding process
The previously described assimilation procedure, performed in situ and in real-time, can be used in the context of welding control. If the stochastic prediction on the quantity of interest T 3|τ * is not satisfying with regards to the criterion Q(T 3|τ * ) = 1, a change in the command u(t) can be implemented as described in "Real-time control" section. This implementation is performed here. In Fig. 15, we show the time evolution of the pdf associated with the prediction on T 3|t , with or without control. In the case without control, the sharp time evolution is due to changes in the pdfs of σ and Pe along the data assimilation steps. We observe that the quantity Q(T 3|τ * ) is much larger than 1, indicating overheating and wasted energy. On the contrary, implementing the control by varying the magnitude u of the heat source enables to reach the criterion Q(T 3|τ * ) = 1 perfectly, and it also speeds up the convergence of the pdf on T 3|t to the target. In Fig. 16, we indicate the evolution of the command variable along the welding process (in terms of corrections δu i at each assimilation time point t i ). We again observe that the feedback loop is effective and quickly (i.e. much before the final time τ * ) leads to an asymptotic regime in which the command remains almost constant (i.e. δu i ≈ 0). We also show in Fig. 16 the map orders which are used along the data assimilation process when the control is performed. This indicates that an order 1 map is still usually sufficient, but that a few more maps with higher order are required compared to the case with no control (where only the first map was order 4). Eventually, we display in Fig. 17 the evolution in time of the overall CPU cost required to implement the feedback loop, which includes both data assimilation and command synthesis steps. As foreseen, this cost is higher during the first assimilation times when the pdfs on parameters σ and Pe significantly evolve (i.e. when much is learnt from measurement data). Once the asymptotic regime is reached in the model updating procedure, the CPU cost is low (< 1 s) which is compatible with real-time contraints for the considered welding application.

Case 2: control of the welding depth with evolving physical process parameters
This second test case has many similarities with the previous one, the control objective still being Q(T 3|τ * ) = 1. Nevertheless, we now take τ * = 100 and we assume that the welding process experiences an unexpected change in the Peclet number value during service (e.g. due to change in the source velocity or material thermal properties), at t = 40. Consequently, the reference parameters values which are now used to get synthetic (noisy) data are: Starting from the same prior distribution of parameters as in the test case 1, sequential data assimilation using Transport Map sampling and PGD is again performed. The minimization problem associated with the computation of the maps is solved with order 1 information on the derivatives, that is a BFGS algorithm with explicit computation of the gradient from the PGD representation. The complexity of the maps (that is the degree of employed Hermite polynomials) is increased until reaching a tolerance of 10 −3 on the variance diagnostic. We represent in Fig. 18 the evolution in time of the marginals on both parameters σ and Pe. Again, we observe that they become thinner with larger maximal pdf values when the number of data assimilation times increases. We also observe that after the change of the reference value for Pe, the data assimilation algorithm is able to detect this change and infers a mean value that slowly tends to the new reference value (even though right after t = 40, the reference parameter value Pe ref = −55 appears in the tail of the pdf). Meanwhile, during this transient regime, it seems that no additional knowledge is brought for the inference of σ as the associated marginals are stagnating. We also show in Fig. 19 the map orders which are used along the data assimilation process. This particularly indicates that an order 1 map remains sufficient to follow the sudden change in the reference value for Pe. From the dynamical updating of model parameters and with respect to the objective, the control of the process with on-the-fly command synthesis is implemented. We show in Fig. 20 the time evolution of the pdf of T 3|t in the case of a controlled welding process. We observe that the control objective is reached even though pdfs of model parameters have not converged yet around the reference parameter values. This illustrates the interest of the control in a stochastic framework, in which uncertainty on the inferred parameters is taken into account in the synthesis of the command in order to make safe decision. We also plot in Fig. 21 the evolution of the command variable u(t) along the process as well as its incremental corrections δu i at each time point t i ; we clearly observe the change in the command when the physical value of the Peclet number drops at t = 40.

Case 3: control of the welding temperature evolution with prescribed time path
In this last test case, the control objective is to make the temperature T 3|t follow a predefined time path, which comes down to imposing the welding history along the process. We   The prescribed evolution curve for T 3|t is shown in Fig. 22 (dashed red line). It is a ramp increase up to t = 20, then a plateau evolution. In our stochastic framework, the command law is designed so that the predicted mean value of T 3|t follows this target evolution. In practice, at each assimilation time point t i , and from the inferred pdfs on model parameters at this time, a command correction δu i is computed so that the prediction on mean(T 3|t i+1 ) coincides with the target value at the next assimilation time point t i+1 . The evolution of T 3|t predicted from the model with reference parameter values, and without any control, is also shown in Fig. 22 (solid black line).
Starting from the same prior distribution of parameters as in the previous test cases, sequential data assimilation using Transport Map sampling and PGD is performed. The minimization problem associated with the computation of the maps is solved with order 1 information on the derivatives, and the complexity of the maps is increased until reaching a tolerance of 10 −3 on the variance diagnostic. We represent in Fig. 23 the evolution in time of the marginals on both parameters σ and Pe. As expected, we observe that they become thinner with larger maximal pdf values tending to reference parameter values along the data assimilation process. The map orders which are used along this process are shown in  Fig. 24; they again indicate that an order 1 is sufficient, except for first assimilation steps where the complexity of the transformation between the reference density and the first posterior densities is higher.
From the dynamical updating of model parameters and with respect to the objective, the control of the process with on-the-fly command synthesis is implemented. We show in Fig. 25 the resulting time evolution of the pdf of T 3|t . We observe that mean(T 3|t ) quite perfectly matches the target evolution. We also plot in Fig. 26 the evolution of the command variable u(t) along the process as well as its incremental corrections δu i at each time point t i . We observe that during the transient phase (ramp evolution of the target), fast modifications in the command are required while command increments tend to zero once the steady-state target regime is reached. Anyhow, this test case shows that the proposed DDDAS strategy is capable of generating complex and effective command laws.

Conclusions
In this work we presented a procedure to build a numerical feedback loop for the control of a fusion welding process from modeling and simulation, while taking uncertainties into account. In order to perform fast computations and permit real-time exchanges between the physical system and its virtual twin, PGD model reduction and Transport Map sam- pling were used in several numerical tasks along the feedback loop. In particular, the explicit dependency on the model parameters inside the PGD model as well as the suitable sampling and integration framework offered by transport maps enabled to effectively perform data assimilation, uncertainty quantification, and predictive control. The implementation of the feedback loop for various control scenarios illustrated the interest and performance of the proposed approach. This approach thus appears to be a relevant tool for real-time feedback control in the DDDAS framework. Future works should focus on the extension of the approach to more complex (e.g. nonlinear) models, associated with modeling errors that may be a priori considered in the Bayesian framework but also a posteriori corrected from data-based learning and enrichment. Dealing with a larger number of model parameters and control variables in the DDDAS context is also a research topic of interest that will be investigated in forthcoming works.