DMC-TPE: tree-structured Parzen estimator-based efficient data assimilation method for phase-field simulation of solid-state sintering

ABSTRACT In phase-field simulations, accurate material parameters are required to quantitatively predict microstructural evolutions. Non-sequential data assimilations enable the estimation of unknown material parameters by minimizing a cost function that represents the misfit between numerical simulation results and time-series observation data. In this study, a new non-sequential data assimilation method Minimizing the Cost function using tree-structured Parzen estimator (TPE), namely DMC-TPE, with higher estimation accuracy than that of a conventional method (DMC-Bayesian optimization (BO)) is developed. The estimation accuracy of DMC-TPE is compared with that of DMC-BO via numerical experiments, where these methods are applied to a phase-field simulation of solid-state sintering. The comparison results demonstrate that the estimation accuracy of DMC-TPE is higher than that of DMC-BO, specifically in cases where numerous parameters have to be estimated, because TPE can continuously minimize the cost function by increasing the number of iterative minimization calculations. Furthermore, DMC-TPE provides less scatter in the estimation results than that in the case of DMC-BO. The DMC-TPE developed herein leads to highly accurate PF simulations of microstructural evolution by simultaneously estimating the states and many unknown material parameters with high accuracies based on experimental data. GRAPHICAL ABSTRACT IMPACT STATEMENT This study developed a novel data assimilation method that can estimate multiple material parameters with high accuracy and less scatter, leading to highly accurate phase-field simulations of solid-state sintering.


Introduction
Phase-field (PF) method [1] is one of the most powerful numerical methodologies for solving free boundary problems.The PF method enables us to simulate complicated evolutions of interfaces without explicitly tracking them.PF simulations are used to predict microstructural evolutions such as during solidification [2][3][4], grain growth [5,6], and sintering.In recent years, various PF simulations of sintering have been performed [7][8][9][10][11][12][13][14] because sintering is used to manufacture various materials including polycrystalline bulk superconducting materials [15,16] and sintered magnets [17,18] and is one of the fundamental technologies for additive manufacturing.To quantitatively predict microstructural evolutions using PF simulation, accurate physical properties and model parameters are required.Physical properties, model parameters, and other material parameters are conventionally obtained via experiments, first-principles calculations, and molecular dynamics simulations.However, the identification of all the material parameters needed for PF simulations using conventional approaches is time-consuming.Additionally, some model parameters are phenomenological, and, therefore, their determination via conventional approaches is challenging.
Data assimilation (DA) is a computational technique that integrates numerical simulation models with experimental data based on Bayesian inference [19].By treating the simulation results and experimental data as stochastic variables that follow a probability density function (PDF), DA can estimate unobservable states and unknown parameters by considering the uncertainties in the simulation models and experiments.Although DA has been mainly developed for meteorology and oceanography [20][21][22], recently, DA has also been employed in material science.In several studies [23][24][25][26][27][28][29][30][31][32], DAs have been applied in PF simulations to estimate experimentally unobservable states and unknown material parameters according to experimental data.
DA methods are generally classified into two types: sequential and non-sequential DAs.Ensemble Kalman filter (EnKF) [33], which is a typical sequential DA method, and EnKF-based DA methods have been utilized in some previous studies [23][24][25][26][27][28][29].EnKF-based DA methods approximate PDFs using numerous simulation results (termed ensemble approximation or Monte Carlo approximation) [33,34].Although the implementation of EnKF-based DA in a PF model is easy, the computational cost is high because a large number of simulations are required for the ensemble approximation.Furthermore, the results of the state and material parameter estimations depend on the accuracy of the ensemble approximation.
Four-dimensional variational method (4DVar) [35,36] is a characteristic non-sequential DA method.In 4DVar, the misfit between the time-series observation data and the simulation results is defined as a cost function.By minimizing the cost function, 4DVar estimates the optimal initial states and parameters (further details are provided in Section 3).4DVar has the advantage of low computational cost.Moreover, various observational data and physical constraints can be used in 4DVar [36,37].Nevertheless, when 4DVar is applied in a PF model, complicated calculations are needed for deriving an adjoint model [36] to obtain the gradient of the cost function, which is required for the variational method.Additionally, if the PF model is even slightly changed, the adjoint model must be recalculated.
To avoid the calculation of the adjoint model, Altaf et al. [38] have proposed a variational DA method that approximates the gradient of the cost function using a reduced-order approach [39,40].This DA method can be easily implemented in numerical simulation models.However, this method still needs to recalculate the gradient approximation when the numerical model is changed.Liu et al. [41,42] have developed an ensemble-based 4DVar (En4DVar) that eliminates the adjoint model by deriving the gradient of the cost function using the ensemble approximation.In a previous study [43], we have applied En4DVar to the PF simulation of solid-state sintering and successfully estimated sintered states and multiple material parameters: two types of atomic diffusion coefficients and two types of mobilities.Nevertheless, because En4DVar utilizes the ensemble approximation, it is computationally expensive and the estimation results depend on the accuracy of the ensemble approximation.
Eliminating the gradient of the cost function from the minimization calculation is one of the simplest approaches to realize a DA method that can overcome the computational cost and implementation problems.In a previous study [44], we have established a nonsequential DA method: Data assimilation method Minimizing a four-dimensional Cost function using Bayesian Optimization (DMC-BO).Bayesian optimization (BO) [45,46] can minimize a black-box function at low computational cost without calculating its gradient.Thus, the application of DMC-BO in numerical simulation models is easy, and no complicated calculations for the adjoint model are required.Furthermore, the states and multiple material parameters can be estimated at low computational costs by applying DMC-BO to a PF model of solid-phase sintering and conducting numerical experiments [44].
As mentioned above, although DMC-BO is an effective DA method, its estimation accuracy still needs to be further improved.In contrast, in the field of mathematical optimization, advanced BO algorithms have been developed to improve performances of various optimization problems.Using a gradient-free superior algorithm gradient for non-sequential DAs, the estimation accuracy of this type of DAs can be improved.
In this study, we developed a novel DA method using the tree-structured Parzen estimator (TPE) to minimize the four-dimensional cost function (DMC-TPE).TPE [47,48] is an optimization algorithm based on the Bayes' theorem and is commonly used for hyperparameter optimization in machine learning.Similar to BO, the TPE minimizes a black-box function without its gradient.Bergstra et al. [47] reported that the TPE decreased an objective function more than that in the case of BO in an optimization problem with 32 variables.Additionally, the TPE has been employed to simultaneously optimize over a hundred parameters [48].Thus, DMC-TPE is expected to lead to more accurate state and parameter estimations than those in the case of DMC-BO, specifically when many material parameters have to be simultaneously estimated, with the same advantages as those of DMC-BO.However, the performances of optimization algorithms depend on the problem to be applied.Therefore, we validated the accuracies and computational costs of the estimations of states and material parameters using DMC-TPE by performing numerical experiments called twin experiments [21,22].Although DMC-TPE can be applied to any PF model, herein, a PF simulation of the solid-state sintering of silver particles is selected as a benchmark problem for validation as in previous studies [43,44].Moreover, we conducted twin experiments using DMC-BO and compared the estimation results of DMC-BO with those of DMC-TPE.The results of these comparisons demonstrate that the DMC-TPE developed in this study is a powerful DA method for estimating material parameters and improving the prediction accuracies of PF simulations.

Phase-field model
In this study, the PF model of solid-state sintering proposed by Wang [7] was applied to DMC-TPE and DMC-BO.This PF model enables the simulation of the microstructural evolution and macroscopic densification of sintered compacts by analyzing the atomic diffusions and rigid-body motions of sintered particles.In this section, we outline the main features of this PF model.Further details of the PF model of solidstate sintering can be found in the literature [7].
To simulate the time evolutions of sintered compacts, we defined two types of PF variables: density field, ρ, and orientation field, η α , (α = 1, 2, 3, . . .N p , where N p is the number of particles).ρ represents the probability of particle existence and is defined as a conserved variable to satisfy mass conservation.η α identifies individual particles and is treated as a non-conserved variable because the volume of each particle varies with grain growth.ρ and η α smoothly changed from 0 to 1 on the surface and in grain boundary regions.
The total free energy F is defined as follows [7]: ( ) dV; ( where f bulk is the chemical free energy density in the bulk. The second and third terms on the right-hand side of Equation ( 1) are the gradient energy densities.κ ρ and κ η denote the gradient energy coefficients that are defined as where γ s is the particle surface energy, γ gb is the grain boundary energy, and δ is the interfacial thickness.The morphological changes of the sintered particles were simulated by calculating the time evolution of ρ.The time evolution equation of ρ was derived by adding an advection term to the Cahn-Hilliard equation [49] to analyze the rigid-body motion of the particle: where v adv,α is the advection velocity of the αth particle and M ρ is the diffusion mobility, which is defined as follows [14]: where M vol , M vap , M surf , and M gb are the volume, vapor, surface, and grain boundary diffusion mobilities of atoms, respectively.These diffusion mobilities are calculated using the diffusion coefficients of the corresponding diffusion paths [14].q(ρ) represents an interpolation function, and q(ρ) = ρ 3 (10-15ρ +6ρ 2 ).Rigid-body motion is assumed to consist of translational and rotational motions.Thus, v adv,k is presented by where v tr,α and v ro,α are the translational and rotational velocities of the αth particle, respectively, which are proportional to the translational and rotational mobilities m tr and m ro , respectively [7].By calculating the time evolution of η α , grain boundary migration was simulated.The time evolution equation of η α was derived by inserting the advection term v adv,α into the Allen -Cahn equation [50]: where M η is the mobility related to grain boundary migration.
In the numerical simulation, we solved the time evolution equations of ρ and η α using the first-order Euler finite difference method (FDM) for time integration and second-order central FDM for space derivatives.

Cost function
We applied DA to the PF model to estimate the changes in the morphologies of sintered particles and grain boundary migration during solid-state sintering from time t = 0 to t = t end .For DA, initially, we defined the augmented state vector x t according to numerical simulation results, where x t contains all state variables in the computational domain at time t (0 ≤ t ≤ t end ) and the material parameters to be estimated.The dimension of x t , denoted as l d , represents the sum of the number of the order variables ρ and η α at all finite-difference grids in the computational domain and the number of parameters to be estimated.Moreover, we define the observation vector y t , which comprised the experimental data.The dimension of y t , represented by m d , indicates the number of coordinates at which the experimental data are acquired.Both x t and y t have numerical calculation and experimental errors.In DA, these errors are assumed to be independent of each other and follow a Gaussian distribution [36].Under these assumptions, x t and y t were treated as stochastic variables that followed a Gaussian PDF.
x t , which was obtained from a numerical simulation, was expressed as where M denotes the model operator.In this study, M corresponds to the PF model of solid-state sintering described in Section 2. The cost function, which represented the magnitude of the data misfit between x t and y t from time t = 0 to t = t end , was derived from the maximum likelihood estimation method [35,36,43]: where x b 0 is an initially estimated augmented state vector, called the background state vector.The superscript T indicates transposition.B and R t are the covariance matrices of the background and observation errors, respectively.B represents the magnitude of the numerical error in the background state, and R t denotes the magnitude of the experimental error in the observational data.In this study, B and R t were set to l d × l d and m d × m d diagonal matrices, respectively, because we assumed that the errors followed Gaussian distributions of mean zero and the cross-correlation between errors was 0. H t is an observation operator that extracts quantities comparable to y t from x t .The optimal initial states and estimated material parameters are acquired as the optimal state vector x a 0 when the cost function J(x 0 ) is minimized.Because the time series x t is necessary to calculate Equation ( 7), a PF simulation must be performed from time t = 0 to t = t end using x 0 as the initial condition for the evaluation of J(x 0 ).
Conventional non-sequential DA methods, including 4DVar, minimize J(x 0 ) using the cost function gradient ∇J(x 0 ) (= ∂J(x 0 )/∂x 0 ) calculated using the adjoint model [36].Nevertheless, the complicated calculation of adjoint models hinders the implementation of 4DVar in the PF model.In contrast, the DMC-TPE proposed in this study and DMC-BO [44] utilize the TPE and BO to minimize J(x 0 ), respectively.Both the TPE and BO do not require ∇J(x 0 ) for the minimization calculation of J(x 0 ).In the next section, we have comprehensively explained the optimization algorithm of the TPE.

Tree-structured Parzen estimator
The TPE is one of the Bayes' theorem-based iterative methods for optimizing a black-box function.Although the TPE is categorized as a type of BO, it uses kernel density estimation, whereas a typical BO employs Gaussian process regression [51].In the DMC-TPE, the function of Equation ( 7) is treated as an objective function for minimization.The objective function outputs y = J(x 0 ) when it receives the augmented state vector x 0 as an input variable.
For N iterative calculations to minimize y using the TPE, we denoted the nth x 0 and y as x 0 (n) and y(n) (n = 1, 2, . . ., N), respectively.Furthermore, we defined the datasets D(n) = {(x 0 (n), y(n))}.Note that N corresponds to the number of data accumulated during minimization.Figure 1 shows a schematic of the minimization of a one-dimensional black-box function using the TPE when N = 10.Hereinafter, the processes labeled (i)-(iv) in Figure 1 are thoroughly described.
The first step in the TPE algorithm is to classify the N datasets (Figure 1(i)) into two sets: L and G (Figure 1(ii)).L contains the data D(n) with y(n) < y*, whereas G comprises the data D(n) with y(n) ≥ y*.y* is a threshold value and is chosen by the hyperparameter τ such that τ = p(y < y*) [47].Practically, L consists of N L datasets selected in ascending order with respect to y [48], whereas G contains (N -N L ) datasets.Although some definitions of N L have been proposed, herein, we used the following definition [52]: The TPE defines conditional PDF called prior distribution, p(x 0 |y), as where l(x 0 ) and g(x 0 ) are the PDFs calculated using L and G, respectively (Figure 1(iii)).In the TPE, x i 0 (1 ≤ i ≤ d, where d is a dimension of the input variable x 0 ) included in the input and x 0 are assumed to be independent of each other.Therefore, l(x 0 ) and g(x 0 ) can be expressed as follows: where l i (x i 0 ) and g i (x i 0 ) are the PDFs calculated presented by the following equations derived from kernel density estimation [53]: where w i (n) j and σ i (n) j (j = L, G) are the weights and bandwidths of the ith input parameter of the nth dataset in the set j, respectively, and are dynamically determined during iterative minimization process [52].K(x| μ, σ) is a kernel function and is defined as Equation ( 14) represents the Gaussian distribution of the mean μ and the variance σ.
To determine the optimal x 0 to minimize y, the TPE utilizes an expected improvement (EI) acquisition function.An input x � 0 that maximizes the EI acquisition function is selected as the candidate x 0 (N +1) to be obtained.The EI acquisition function is expressed by the following equation [47]: where p(y |x 0 ) is the posterior distribution obtained from the Bayes' theorem: Here, using τ, p(x 0 ) can be expressed as By substituting Equations ( 9), (16), and (17) into Equation ( 15), we obtained the following acquisition function: Equation (18) demonstrates that the maximization of l (x 0 )/g(x 0 ) maximizes the EI acquisition function.Therefore, the TPE algorithm yields the candidate x 0 (N +1) according to the evaluation of l(x 0 )/g(x 0 ) (Figure 1(iv)).

Flowchart of DMC-TPE
The flowchart for the estimation of the optimal augmented state vector x a 0 using DMC-TPE is described as follows (Figure 2): (iii) calculations of the probability density functions l(x 0 ) and g(x 0 ) using kernel density estimation; and (iv) determination of the candidate x 0 (N +1) according to the maximum of l(x 0 )/g(x 0 ).
Step 1: Set the range of possible values for the state variables and material parameters included in x 0 .
Step 2: Determine x b 0 using arbitrary initially estimated state variables and initially estimated parameters.
Step 3: Create the initial datasets D(n) = {x 0 (n), J(x 0 (n))} (1 ≤ n ≤ N).Inputs of the datasets are the arbitrary augmented state vectors x 0 (n), which are defined by the state variables and material parameters within the domain set in Step 1.The outputs of the datasets J(x 0 (n)) are obtained by performing PF simulations of solid-state sintering and calculating the cost function using Equation (7).
Step 4: Classify the datasets into two sets, L and G, based on Equation (8).
Step 7: Perform a PF simulation of solid-state sintering with x 0 (N +1) as the initial condition and then calculate J(x 0 (N +1)) according to Equation (7).
Step 9: Iterate Steps 4-8 while updating N until the end condition is satisfied.The optimal estimate x a 0 is acquired as x min 0 , which affords the minimum cost function obtained during the iterative calculations.
Steps 4, 5, and 6 correspond to the TPE procedures (ii), (iii), and (iv) shown in Figure 1.To construct a new dataset, a PF simulation of solid-state sintering needs to be performed because the cost function must be evaluated.Thus, in DMC-TPE, N is equal to the number of PF simulations required to estimate the states and material parameters.

Concept of twin experiments
To validate the estimation accuracy of the developed DMC-TPE, we conducted numerical experiments called twin experiments using the PF simulation of solid-phase sintering as a benchmark problem.The estimation accuracy of DMC-TPE was validated by comparing the results of the twin experiments achieved using DMC-TPE with those obtained using DMC-BO [44].Silver was selected as the sintered material because its physical properties are relatively isotropic and some physical values are known [54][55][56][57][58].In the twin experiments, at first, we performed a PF simulation of solid-state sintering using an a priori assumed-true initial state and material parameters to create synthetic observation data.Assuming that only morphological changes could be experimentally observed, only the time-series three-dimensional distributions of ρ were used for the synthetic observation data.Additionally, the initial states of the particle shape and grain boundary distributions were assumed to be known because they could be determined [59,60] or deduced [61] from the experimental data.By comparing the a priori assumed-true states and material parameters with the estimated states and material parameters, we investigated the estimation accuracies of DMC-TPE and DMC-BO.The computational cost was evaluated by referring to the number of iterations N until a small cost function was acquired.Furthermore, we examined the dependences of the estimation results on the randomly created initial dataset and random numbers used in optimization (hereinafter, the random elements are collectively referred to as 'randomness') by performing five runs of the twin experiments for each method using different random seeds and evaluating the standard deviations (SDs) of the cost function and estimated parameters.
In the twin experiments described in Section 5.1, four material parameters, that is, diffusion coefficient of silver atoms on particle surfaces, D surf , diffusion coefficient of silver atoms at grain boundaries, D gb , translational mobility, m tr , and mobility related to grain boundary migration, Μ η , were simultaneously estimated, as in previous studies [43,44].These parameters strongly affect the microstructural evolutions of sintered compacts: m tr and Μ η determine the rate of densification and grain growth, respectively; D surf and D gb are typically larger than D vol and D vap .In particular, the PF simulation of solid-state sintering is more sensitive to D gb and m tr , than D surf and M h .As stated in Section 5.2, we increased the number of parameters to be estimated to further validate the estimation accuracy of DMC-TPE.We assumed that D surf and D gb were presented by the Arrhenius law and estimated the pre-exponentials D 0 surf and D 0 gb and the activation energies Q surf and Q gb .Additionally, the grain boundary energy γ gb was evaluated.Thus, in Section 5.2, we estimated seven material parameters: D 0 surf , Q surf , D 0 gb , Q gb , m tr , Μ η , and γ gb .

Conditions for the PF simulation of solid-state sintering
Figure 3 shows the initial distributions of the particles and grain boundaries used in the twin experiments.The computational domain was a cube with sides of 12.8 μm.The spacing of the finite-difference grid was Δx = Δy = Δz = 0.2 μm, and δ = 3Δx.Herein, 12 spherical silver particles were placed in the computational domain as an initial state.The grain boundary regions are described by the quantityΨ, which is defined as The regions with Ψ = 0.5 correspond to the grain boundary.The zero Neumann condition was applied to the boundaries of the computational domain.We calculated the time evolutions during the solid-state sintering of the silver particles for 40 s.
The sintering temperature was assumed to be constant at 600 °C in Section 5.1, while the temperature was assumed to increase at a constant rate of 10°C/s from 300°C to 700°C in Section 5.2 because the preexponentials and activation energies of the diffusion coefficients were simultaneously estimated.Surface energy, γ s and molar volume, V m , were set as γ s = 1.21 J m −2 [54] and V m = 10.27 × 10 −6 m 3 mol −1 , respectively.Volume diffusion coefficient, D vol was set as D vol = 3.00 × 10 −16 m 2 s −1 for Section 5.1 and it was determined by the Arrhenius law with the preexponential D 0 vol = 6.70 × 10 −5 m 2 /s and the activation energy Q vol = 1.90 × 10 5 J/mol for Section 5.2 [56].The temperature dependences of the other parameters used in Section 5.2 were neglected because we focused on investigating the estimation accuracies of DMC-TPE when the number of parameters to be estimated is increased.Diffusion coefficient of silver in the vapor phase, D val , was fixed at 5 × 10 −17 m 2 s −1 [44], which was sufficiently smaller than the surface diffusion or grain boundary diffusion coefficients because the sintering temperature was considerably lower than the melting point of silver.

Definitions for data assimilation
State and observation vectors used in the twin experiments are described hereinafter.In this section, the vectors used for estimating four material parameters, as described in Section 5.1, are defined.The concept of defining the vectors remains the same regardless of the parameters to be estimated.x t is expressed as follows: where ρ t and η t are vectors containing the state variables ρ and η α at all finite-difference grid points at time t, respectively.y t is defined as In the twin experiments, we assumed that the observation data were acquired from the same positions as the finite-difference grid points of the PF simulation.Therefore, ρʹ t is a vector consisting of ρ at all finitedifference grid points for synthetic observation data.Furthermore, to calculate the differences between x t and y t using the same physical quantity, the observation operator H t is defined as follows: Because the initial states of the particle shape and grain boundary distributions were assumed to be known, the initial state vector x 0 (n) optimized by the TPE contained only the material parameters to be estimated: D surf (n), D gb (n), m tr (n), and M η (n),

Computational environment
Implementation of DMC-BO in the numerical model is easy because not only DMC-BO avoids the calculation of the cost function gradient but also various computational libraries can be used for BO (e.g.GPyOpt [62]).For DMC-TPE, a Python library called Optuna [52] can be employed to facilely run the TPE.
In this study, we utilized Optuna (version 2.10.0) to perform the TPE.Thus, DMC-TPE could be applied in numerical models as easily as DMC-BO.In Optuna, the EI acquisition function is maximized by evaluating 24 candidates of l(x 0 )/g(x 0 ), selected based on l(x 0 ) [52].This simplified method leads to good optimization [63,64] while reducing the computational cost for the maximization calculation of l(x 0 )/g(x 0 ) in the highdimensional space.
The DMC-TPE algorithm runs in Python, whereas the PF simulations of solid-state sintering operate in Fortran, which is superior in terms of computational speed.For the PF simulations, we performed parallel computations using four central processing units (Intel Xeon Gold 6150) and four graphics processing units (NVIDIA Tesla V100) with the message passing interface library [65] and a compute unified device architecture library [66].The time required to perform a single iterative minimization was approximately 50 min.

Estimations of states and four material parameters
We conducted twin experiments to estimate four material parameters (namely, D surf , D gb , m tr , and Μ η ) using DMC-TPE and DMC-BO.In this study, for both DMC-TPE and DMC-BO, randomly created initial datasets and random numbers were used.Therefore, the estimation results achieved using both methods depended on randomness.As described in Section 4.1, we validated the dependences of the estimated results on randomness by investigating the SDs of the cost function and the four material parameters to be estimated obtained via five runs of the twin experiments.Table 1 presents the a priori assumedtrue values for the four material parameters (D true surf , D true gb , m true tr , and M true η ).D true surf and D true gb were determined from the literature [57,58].The values of m true tr and M true η were assumed because their actual values were unknown.The assumed values were set such that the morphology of the sintered particle changed with rigid-body motion and grain boundary migration during the PF simulations under the simulation conditions described in Section 4.2.The PF simulations were performed using γ gb = 0.79 J m −2 [55].
Figure 4 shows the PF simulation results acquired using the a priori assumed-true values provided in Table 1.These results demonstrated that neck growth, pore annihilation, and grain boundary migration were simulated by the PF simulation of solid-state sintering for 40 s.In the twin experiments, only the threedimensional time evolution of ρ shown in Figure 4(a) was employed as the synthetic observation data.The time interval of the synthetic observation data was fixed at 2 s, which enabled the estimation of the states and material parameters with high accuracies [43].Thus, 20 synthetic observation data were Table 1.A priori assumed-true values of the four material parameters to be estimated.used to estimate the states and material parameters using DMC-TPE and DMC-BO.Table 2 presents the means and SDs that define the Gaussian distributions of the initially estimated material parameters.The mean values were set to approximately half of the a priori assumed-true values [43,44].The SDs were also set to the same values as those tuned for the highly accurate estimations of the states and material parameters [43].
To demonstrate the difference between the results of the PF simulations obtained without DA and the synthetic observation data, the PF simulation results were obtained using the means of the initially estimated parameters (Figure 5).The particles shown in  1 under isothermal conditions (600°C).(a) the region of ρ > 0.5 is shown as the particle shape.(b) Distributions of the grain boundaries on the cross-section of y = z are described using Ψ, calculated by Equation (19).The evolution of ρ shown in (a) was employed as the synthetic observation data for twin experiments.Because we assumed that the initial particle shape and grain boundary distribution were known, B was defined as B = diag(0, 0, . .., 0, (σ(D surf )) 2 , (σ(D gb )) 2 , (σ(m tr )) 2 , (σ (M η )) 2 ), where σ(•) are the SDs provided in Table 2.Moreover, we assumed that the errors included in the observation data were uniformly σ y at all times, and R t = σ y 2 I, where I is an identity matrix with dimensions m d × m d .In this study, we set σ y = 0.1 [44].
In DMC-BO, the cost function shown in Equation ( 7) is minimized using the BO algorithm.Steps 4-6 of the flowchart shown in Section 3.3 are replaced by performing Gaussian process regression, calculating an acquisition function, and maximizing the acquisition function, respectively.See Ref [44] for further details of DMC-BO.Because the TPE algorithm is based on the EI function as shown in Equation ( 18), the BO algorithm also used the EI acquisition function in this study.GPyOpt maximizes the EI function with the limited-memory Broyden-Fletcher-Goldfarb-Shanno algorithm for bound constrained optimization (L-BFGS-B) method [62,67].There is no difference in the PF model of solid-state sintering used for DMC-TPE and DMC-BO.
Figure 6 shows the transitions of the mean and SD of the minimum of the cost function J(x min 0 ) for N. J (x min 0 ) for N ≤ 5 was evaluated from the randomly created initial dataset, whereas that for N ≥ 6 was obtained by minimization iterations using DMC-TPE and DMC-BO.Hereinafter, J(x min 0 ) values acquired using DMC-TPE and DMC-BO are referred to as J TPE and J BO , respectively.
DMC-BO substantially decreased the mean of J BO by N = 15, and then, J BO was updated less frequently.In contrast, the mean of J TPE was larger than that of J BO for N = 15, whereas J TPE continuously decreased for N > 15.Therefore, the mean of J TPE for N = 21 decreased to a value comparable to that of J BO for N = 15.Furthermore, for N ≥35, the mean of J TPE was consistently smaller than that of J BO .Figure 7 shows a comparison between the (a) initially estimated states, (b) estimated states for N = 15 acquired using DMC-BO, and (c) estimated states for N = 21 obtained using DMC-TPE at t = 40 s.The estimated states shown in Figure 7(b,c) were obtained from one of the results of the five twin experiments in which J TPE or J BO was closest to their respective mean.ε err values in the cases of Figure 7(b,c) were lower than those in the case of Figure 7(a).The grain boundary distributions shown in Figures 7(b,c) reproduced the synthetic observation data shown in Figure 4.These results demonstrated that DMC-TPE estimated the states with accuracies comparable to those of DMC-BO without significantly increasing the computational cost as compared to that for DMC-BO.Table 3 provides the mean and SDs of the estimated material parameters obtained using DMC-BO for N = 15 and DMC-TPE for N = 21.Both methods accurately estimated D gb and m tr with small SDs among the four material parameters because D gb and m tr are sensitive parameters in terms of the time evolutions of the state than those of D surf and M η in the PF simulation of isothermal sintering [43]; that is, accurately estimated D gb and m tr led to accurate estimation of the sintered states.Although a previous study [44] has reported that DMC-BO affords highly accurate estimation results with less scatter for D gb and m tr , this study indicates that DMC-TPE also estimates the sensitive parameters with high accuracies and small SDs.
As mentioned above, DMC-TPE more frequently updates the minimum of the cost function as compared to the case of DMC-BO.Hereinafter, we discuss the differences between the iterative calculations of minimizing the cost function conducted using DMC-TPE and DMC-BO.Figure 8 shows histograms of the frequencies of the candidate values for each material parameter suggested by DMC-BO and DMC-TPE.The histograms were calculated using the results of the five twin experiments.For DMC-BO, the candidate values are frequently presented at the upper and lower bounds of the search range, specifically for N ≥ 20.This 'boundary search' occurs because the BO algorithm performs global searches by referring to the regions with less data.Note that the boundary search takes place not only in DMC-BO but also in other numerical algorithms using BO (e.g.Ref [68]).J BO rarely decreases during the boundary search when the true value cannot be found in the boundary region, for instance, in the cases of the present twin experiment and actual situation where the ranges are set sufficiently wide to include unknown true values.Therefore, DMC-BO less frequently updates J BO for N ≥ 20.Furthermore, the SD of J BO is unlikely to decrease during a global search because a global search is more sensitive to randomness than local optimization, which determines candidate values based on the best data.The global search performed in DMC-BO can be suppressed by setting a hyperparameter to encourage local optimization; nevertheless, in these cases, the candidates can easily fall into a local optimum solution.
In contrast, for DMC-TPE, all histograms shown in Figure 8(b) exhibit unimodal distributions with peaks  close to the true values.The unimodal distributions indicate that the TPE performed local optimizations with parameters that are closer to the true optimal solution than the parameters suggested by BO.As described in Section 3.2, DMC-TPE determines the next candidate by referring to the set l(x 0 ) consisting of data with small cost functions.Additionally, l(x 0 ) is evaluated under the assumption that the parameters to be optimized are independent of each other (Equation ( 12)).This assumption facilitates the consideration of the effect of each material parameter on the decrease in the cost function.Therefore, a value in the boundary region is rarely suggested, and the cost function in the case of DMC-TPE continuously decreases with an increase in N. The SD of J TPE also decreases because the candidates gradually converge to the true value.It should be noted that the evaluation of the effect of each parameter on the cost function in DMC-TPE led to the larger SDs for the nonsensitive parameters (D surf and M η ) obtained using DMC-TPE than those obtained using DMC-BO.Moreover, TPE prevents the falling of candidates into the local solution by considering not only l(x 0 ) but also g(x 0 ).These results demonstrate that DMC-TPE is superior to DMC-BO as a gradientfree non-sequential DA method.DMC-TPE can estimate the states and material parameters with higher accuracies and smaller dependences on randomness with an increase in N, than those in the case of DMC-BO.
In the four parameter estimations, the significant differences in estimation accuracy between DMC-TPE and DMC-BO were not clearly shown because J TPE and J BO were reduced with a relatively small number of iterations.Thus, in the next section, the twin experiment for the seven parameter estimations was conducted, which is expected to require many iterations to reduce the cost function.

Estimations of states and seven material parameters
As described in Section 1, the TPE is effective for problems in which numerous parameters have to be simultaneously optimized.Therefore, in this section, we simultaneously estimated the state and seven material parameters using DMC-TPE and DMC-BO applied to the PF simulations of non-isothermal solid-state sintering, in which the temperature was increased from 300°C to 700°C.As explained in Section 4.1, the following seven parameters were assessed: were not different from those provided in Table 1.
The following ranges of possible values for D 0 surf , Q surf , D 0 gb , Q gb , and γ gb during iterative calculations for n times were determined: 0 ≤ D 0 surf (n) ≤ 2.0 × 10 2 m 2 s −1 , 2 × 10 5 ≤ Q surf (n) ≤ 2.5 × 10 5 J mol −1 , 0 ≤ D 0 gb (n) ≤ 2.5 × 10 −5 m 2 s −1 , 8.0 × 10 4 ≤ Q gb (n) ≤ 1.0 × 10 5 J mol −1 , and 0.5 ≤ γ gb (n) ≤ 1.2 J m −2 .The maximum possible values of D 0 surf and D 0 gb were approximately twice the corresponding true values.The ranges for Q surf and Q gb were determined by trial and error such that numerical divergence did not occur.The range of γ gb was determined under the assumption that γ gb < γ s .Table 4 provides the means and SDs used to define the Gaussian distributions of the seven initially estimated material parameters.To set small initially estimated values for D surf and D gb , the initially estimated values of D 0 surf and D 0 gb were fixed at approximately half of their true values, whereas those of Q surf and Q gb were set to the maximum values in their respective ranges.Although tuning the SDs is an important factor for achieving highly accurate estimation results [39], the SDs used herein are not elaborately tuned because we aim to investigate the difference between the estimation results of DMC-TPE and DMC-BO.Other calculation conditions are the same as those stated in Section 5.1.
Figure 9 shows the transitions of the means and SDs of J TPE and J BO .The mean of J TPE decreased with an increase in N and was smaller than that of J BO for N ≥8.The differences between the means of J TPE and J BO for N = 100 were larger than those obtained from the estimation results of four parameters.As described in Section 5.1, DMC-TPE can avoid the boundary search.The larger the number of parameters to be estimated (i.e. the dimension of the search space), the stronger the influence of the frequency of boundary search on the minimization performance degradation.Therefore, in the simultaneous optimization of several material parameters, DMC-TPE works more effectively than DMC-BO.Figure 10(a) shows the initially estimated state obtained at t = 40 s (700°C).Figure 10(b,c) show the estimated states for N = 100 obtained from DMC-BO and DMC-TPE, respectively.To clearly indicate the difference between ε err for DMC-BO (b) and DMC-TPE (c), the distributions of ε err in the computational domain are shown in the second row of Figure 10.Although both DA methods reduced ε err , the state estimation result achieved using DMC-BO left regions with larger ε err than that of DMC-TPE.This reflects the difference between J TPE and J BO for N = 100, as shown in Figure 9.These results demonstrated that DMC-TPE estimated the states with higher accuracies than those of DMC-BO with an increase in the number of parameters to be estimated.Moreover, the SD of J TPE was less than half that of J BO for N ≥19, and for N = 100, the SDs of J TPE and J BO differed by one order of magnitude (1.4 × 10 3 and 1.1 × 10 4 , respectively).Therefore, DMC-TPE provides stable and accurate estimation results of states with less dependence on randomness.
Table 5 presents the material parameters estimated using DMC-TPE and DMC-BO for N = 100.DMC-TPE estimated D 0 surf , Q surf , and m tr with higher Means of the initially estimated parameters 0.50 × 10 2 2.50 × 10 5 6.00 × 10 −6 1.00 × 10 5 5.00 × 10 −11 1.00 × 10 −7 1.00 SDs of the errors included in the initially estimated parameters 0.25 × 10 2 5.00 × 10 3 3.00 × 10 −6 1.00 × 10 3 5.00 × 10 −12 2.00 × 10 −8 0.20 accuracies and smaller SDs than those of DMC-BO.The SDs of m tr and g gb , obtained using DMC-TPE, were much smaller than those of DMC-BO.For D 0 gb and Q gb , which yield the sensitive parameter D gb , there was no significant difference in the estimation accuracy between DMC-TPE and DMC-BO, with DMC-TPE being slightly more accurate and less scatter.On the other hand, the SD of M h obtained using DMC-TPE was larger than that of DMC-BO.As mentioned in the previous section, this result reflects the characteristics of DMC-TPE algorithm and indicates that M h is a non-sensitive parameter.From these results, DMC-TPE demonstrated that it provides accurate and less scattered estimation results for sensitive parameters without the calculation of the cost function gradient, even with an increase in the number of parameters to be estimated.

Conclusions
This paper proposed a novel DA method called DMC-TPE that utilized the TPE algorithm for the minimization calculation of the cost function.DMC-TPE can be easily implemented in PF models because it does not require the gradient of the cost function.In this study, we performed twin experiments where DMC-TPE and DMC-BO were applied to the PF simulations of solidstate sintering to investigate the differences between the estimation accuracies and computational costs of these two methods.The following details were obtained from the twin experiments: (1) When four material parameters were simultaneously estimated, both DMC-TPE and DMC-BO estimated the states with high accuracies at  γ gb [J m −2  parameters with high accuracies, a large N was required.In these cases, DMC-TPE estimated the states and parameters more accurately than DMC-BO because it continuously decreased the cost function with an increase in N.
(3) DMC-TPE decreased the SD of the cost function with an increase in N.For N = 100 in the seven parameter estimations, the SD obtained from DMC-TPE was approximately one-tenth of that of DMC-BO.These results indicated that DMC-TPE estimated the states and parameters with less dependence on the randomly created initial dataset and random numbers used in optimization.
From the above results of twin experiments, the advantages of DMC-TPE are summarized as follows: • Implementation of DMC-TPE in simulation models is easy.• DMC-TPE continuously reduces the cost function with an increase in N. • The decrease in the cost function is less dependent on the randomness.• The estimation accuracy of DMC-TPE is more accurate than that of DMC-BO, particularly when the number of parameters to be estimated is large.
Based on the results obtained herein, we conclude that DMC-TPE is a powerful DA method for the simultaneous estimation of unobservable states and many material parameters.While the advantage of DMC-TPE was demonstrated, further validations are required to use actual experimental data for DMC-TPE (e.g. the sufficient number of observation data for accurate estimation when the time interval of observation data is long, the determination method of R t when experimental data have non-uniform errors, and the selection of which parameters to make temperature dependent when experiments are conducted under non-isothermal condition).Therefore, we plan to conduct further validation of DMC-TPE in future studies.The clarification of unknown material parameters related to solid-state sintering by integrating in situ experimental observation data and PF simulation results of solid-state sintering using DMC-TPE will have a significant impact on the prediction and elucidation of microstructural evolution and material parameters, leading to further advances in material science.

Disclosure statement
No potential conflict of interest was reported by the author(s).

Figure 1 .
Figure 1.Schematic of the minimization of a one-dimensional black-box function using TPE: (i) initial datasets; (ii) classification of the initial dataset into two sets: L and G;(iii) calculations of the probability density functions l(x 0 ) and g(x 0 ) using kernel density estimation; and (iv) determination of the candidate x 0 (N +1) according to the maximum of l(x 0 )/g(x 0 ).

Figure 2 .
Figure 2. Flowchart for estimating an optimal state vector, x a 0 , using DMC-TPE.

Figure 3 .
Figure 3. Initial distributions of particles and grain boundaries used for the PF simulations of solid-state sintering of silver particles.(a) The region of ρ > 0.5 is shown as the particle shape.(b) Distributions of grain boundaries on the cross-section of y = z are described using Ψ calculated by Equation (19).

Figure 4 .
Figure 4. Snapshots of state evolutions obtained by the PF simulation using the a priori assumed-true values provided in Table1under isothermal conditions (600°C).(a) the region of ρ > 0.5 is shown as the particle shape.(b) Distributions of the grain boundaries on the cross-section of y = z are described using Ψ, calculated by Equation(19).The evolution of ρ shown in (a) was employed as the synthetic observation data for twin experiments.

Table 2 .Figure 5 .
Figure 5. Snapshots of state evolution obtained via PF simulation using the initially estimated parameters presented in Table2.(a) the region of ρ > 0.5 is shown as the particle shape that is colored according to the absolute error ε err .(b) Distributions of the grain boundaries on the cross-section of y = z are described using Ψ, calculated by Equation(19).

Figure 5 (
Figure5(a) are colored based on an absolute error, ε err , which was calculated by ρ for the simulation results and synthetic observation data at each finite-difference grid point.The PF simulation results acquired without DA were significantly different from the synthetic observation data, particularly on the neck and where the pores were annihilated.Because we assumed that the initial particle shape and grain boundary distribution were known, B was defined as B = diag(0, 0, . .., 0, (σ(D surf )) 2 , (σ(D gb )) 2 , (σ(m tr )) 2 , (σ (M η ))2 ), where σ(•) are the SDs provided in Table2.Moreover, we assumed that the errors included in the observation data were uniformly σ y at all times, and R t = σ y

Figure 6 .
Figure 6.Transitions of the means and SDs of the minimum of the cost function J(x min 0 ) obtained from the estimation of the four parameters using DMC-BO and DMC-TPE.The means and SDs were calculated based on the results of five twin experiments for each method.

Figure 7 .Table 3 .
Figure 7.Comparison between the (a) initially estimated states, (b) estimated states for N = 15 obtained using DMC-BO, and (c) estimated states for N = 22 acquired using DMC-TPE at t = 40.The upper row shows the region of ρ > 0.5 as the particle shape that is colored according to the error ε err .The lower row shows the grain boundary distributions on the cross-section of y = z described using Ψ, calculated by Equation(19).

Figure 8 .
Figure 8. Histograms of the frequencies of the candidate values for each material parameter suggested by (a) DMC-BO and (b) DMC-TPE.Dashed lines indicate the true values for each material parameter.

Figure 9 .
Figure 9. Transitions of the means and SDs of the minimum of the cost function J(x min 0 ) obtained from the estimation of seven parameters using DMC-BO and DMC-TPE.The means and SDs were calculated via five twin experiments for each method.

Figure 10 .
Figure 10.Comparison between the (a) initially estimated states obtained by PF simulation using the initially estimated parameters presented in Table4and the estimated states at N = 100 acquired using (b) DMC-BO and (c) DMC-TPE at t = 40 s (700 °C).The first row demonstrates the regions of ρ > 0.5 as the particle shape that is colored according to the error ε err .The second row shows the distributions of ε err .The third row shows the grain boundary distributions on the cross-section of y = z described using Ψ, calculated by Equation(19).

Table 5 .
Means and SDs of the seven parameters estimated using DMC-BO and DMC-TPE.

Table 4 .
Means and SDs of the Gaussian distributions used to determine the seven initially estimated parameters.
low computational costs.For the parameter estimation, D gb and m tr , which strongly affected the state evolution, were particularly estimated with high accuracies and less scatters.(2)To simultaneously estimate states and seven