MIOpt: optimization framework for backward problems on the basis of the concept of materials integration

ABSTRACT In materials design, it is very difficult to accurately design forward problems owing to the variety of scales and phenomena to be considered and the increasing number of input and output variables. In addition, it is necessary to reduce the number of computationally expensive calculations to perform backward problem to obtain the optimal set of parameters for the input to a forward problem. In this paper, we describe the development of the MIOpt (MI Optimizer), a framework for backward problem. The process – structure – property – performance relationship (PSPP relation) of materials is modeled as forward problems in MInt, and constructed in terms of materials engineering concepts. This can be linked to backward analysis via MIOpt, which greatly improves the efficiency of trial – and – error materials design on a computer. In addition, two proprietary optimization algorithms are implemented, allowing the use of an extension of Bayesian optimization and an update algorithm for machine learning models, which reduces the actual time required for optimization. The benchmark problems were the minimum search problem for the analytical function and the heat input condition optimization problem for creep rupture time in welds. Efficient Bayesian optimization method was confirmed to be faster and more stable than the conventional sequential Bayesian optimization method. GRAPHICAL ABSTRACT IMPACT STATEMENT In this study, we would like to report on the development of MIOpt and j-parallel Bayesian Optimization (BO) algorithm. MIOpt is an optimization framework of forward problem workflow design system (MInt) for structural materials and connected directly to MInt, that enables us to calculate backward problem more easily. The j-parallel BO algorithm enables optimization more robust to get convergence.


Framework for forward problem design and backward problem of structural materials
With the recent increase in the number of explanatory variables and the diversification of objective variables in materials design, the search space has increased, making it difficult to efficiently utilize conventional design methods.On the other hand, with the advancement of computational techniques such as first-principles calculations, thermodynamic equilibrium calculations, phase field methods, and finite element methods as elemental technologies, the performance of structural materials can be considered effective according to the process -structure-property -performance (PSPP) models [1].This is being actively discussed in movements such as integrated computational materials science (ICME) [2] and materials genome initiative (MGI) [1].In Japan, the idea of materials integration has been proposed [3,4] and realized as a computer system materials integration by network technology (MInt) that combines experiment, theory, computational science and data science [5,6].
In materials development, a backward problem is required to design the optimal process on the basis of the intended performance; although the P-S-P-P direction is the forward direction as a causal law, it is highly nonlinear, and it is necessary to assume that the response surface is multimodal.There are several commercial packages for optimization (Isight [7], and HEEDS [8]), however, there is no environment in which we can perform backward problem that is directly linked to the forward problem workflow in materials engineering.To meet these requirements, we have developed an optimization platform of Materials Integration Optimizer (MIOpt).MIOpt is a computational infrastructure that links with MInt, specifies explanatory variables and objective functions expressed in materials descriptors, and can be connected to backward problems.In addition, two algorithms were implemented and validated, which improved the efficiency of Bayesian optimization (BO) and update machine learning models.

Linking to MInt and MIOpt
MInt is a system for solving structural material problems by designing a workflow as shown in Figure 1, expressing the relationship between prediction models and input/output.MInt consists of three components: (1) a workflow design and execution function that describes the relationship between prediction models, (2) meta-information and vocabulary management system for prediction models, and (3) experimental and calculation result database.
Considering the use of input and output data as training data for parameter surveys and machine learning model developing, users are encouraged to organize information in accordance with materials engineering information.Forward problems designed in this way can be interconnected by controlling the connection ports.As a result, a variety of workflows can be constructed.
Basis of the PSPP concept, a number of workflows have been developed to consistently predict the performance of structural materials [6,[9][10][11][12][13][14][15][16][17][18].In addition, application programing interface (API) to utilize MInt functions from external programs were developed as MInt-API, such as controlling workflow execution, registering and updating meta-data, and retrieving data from databases.Security measures are taken by authenticating with access tokens.
In this study, we have utilized these APIs to build a framework for backward problem, MIOpt, which can utilize information from the forward problems handled by MInt to optimize the input parameter set for structural materials forward problems shown in Figure 2. In MIOpt, a user selects a workflow to be optimized linking to MInt and sets the explanatory and objective variables from the input -output descriptors in the workflow.Furthermore, from the configuration screen of MIOpt shown in Figure 3, the search range is set, and the optimization calculation is executed utilizing algorithms of mathematical optimization.The optimization algorithms are implemented with common algorithms such as linear regression, the steepest descent, and BO.In addition, we have implemented two algorithms in this study: an algorithm extended to improve the efficiency of BO while making effective use of the computational environment and an algorithm for updating the machine learning model.Details are described in Sections 3.2, 3.3, respectively.

Optimization test
Here, the Rosenbrock function, an analytic function, is taken as the benchmark function.This function having two variables, x 1 and x 2 , is defined by Equation (1) with a global minimum of 0 at ðx 1 ; x 2 Þ ¼ ð0; 0Þ: The functional form shown in Figure 4(a) is simple, but it must be well controlled near the variable that takes the minimum value because the region around the coordinate that takes the minimum value is smooth.First, we designed a workflow for the forward problem of the Rosenbrock function.Then, we performed BO using the GP -EI (Gaussian Process -Expected Improvement) method, and the results shown in Figure 4(b) indicate that the optimization converged to the optimal solution of 0.

Issues and suggestions for improvement of efficiency of BO
In recent trends of utilizing data in materials engineering, known as in materials informatics, BO is frequently used for black-box function optimization and hyperparameter tuning in machine learning model building, as useful tools COMBO [19], PHYSBO [20].In BO, as shown in Figure 5, next calculation points are generated by a Gaussian process and an acquisition function, expressed with the expected value μ and standard deviation σ.Many types of acquisition function are proposed, PI (probability of improvement), EI (expected improvement), and UCB (upper confidence bounds).Emphasizing exploration leads to a greater number of calculations to obtain the optimal solution, whereas emphasizing utilization leads to the possibility of being trapped in a local solution.Therefore, for many acquisition functions, it is necessary to determine a hyperparameter to maintain the balance between exploration and exploitation.Much effort has been made to determine a hyperparameter.The �-greedy method in reinforcement learning involves random search and determining the balance from the value of the action (Q) function [21].In environments with sufficient available computer  resources that are efficient, the optimization methods include parallel Thompson sampling [22,23], batch BO [24,25], evaluation using Expected Improvement (EI) [26,27], and Gaussian process -UCB acquisition function selection methods [28] were proposed.Many of which have been considered for problems involving exploration and exploitation tradeoffs.Although such parallel BO algorithms exist, it is difficult to achieve a balance between exploration and exploitation, which is a major challenge for improving the efficiency of the computation.
In addition, the commonly used sequential BO repeats the evaluation of input parameters for the next calculation on the basis of the results of the previous calculation, so even if there are sufficient computer resources, they cannot all be used, and as a result, optimization may take a long time.
There are also some important BO frameworks developed for materials engineering.CAMEO [29] attempts to predict the next optimal point to evaluate Ge-Sb-Te phase-change memory materials and to expand the knowledge of the system (explore the design space) before optimization.BAREFOOT [30] offers pure BO, and MultiFidelity (sequential BO), or a combination of both.Dragonfly [23] is a library to ensure robustness for global search.The balance between exploration and exploitation is often a tuning parameter in these libraries as well, and if this tuning parameter is not appropriate, the optimization will take a long time.Therefore, in this study, we implemented an algorithm that improves the efficiency of optimization by executing multiple tuning parameters in BO in parallel, and we implemented an algorithm to update a machine learning model.

Execution of BO while balancing between exploration and exploitation simultaneously
The hyperparameters related to the balance between exploration and exploitation are given different names depending on the acquisition function, and in the case of EI, exploration and exploitation are controlled by the hyperparameter jitter (hereafter, j).In GPyOpt [31], which we used in this study for BO, EI is represented by Equation (2).
where f ðx À Þ denotes the best evaluation result recorded until evaluation i (the current best).EI j denotes the expected improvement of best performance if model x is observed, and E is the expectation of posterior distribution of f ðxÞ.The larger the j, the more exploratory the acquisition function becomes.
Therefore, in this study, we investigated an exploration and exploitation-oriented calculation method by using multiple j's simultaneously, and we synthesized the calculation results after each iteration.The pseudo-code for the algorithm of BO using j in parallel (hereafter, j-parallel BO) in this study is shown in Figure 6.First, multiple j's are selected and multiple calculations of the objective function are performed simultaneously for the input parameters presented the  basis on the Bayesian inference.Once the calculations for each j are completed, all the data are synthesized and used as initial data.This procedure is repeated a predetermined number of times.

Optimization by updating the Cubist model
In addition to the algorithm for j-parallel BO proposed in Section 3.2, we constructed an algorithm for updating a machine learning model.The pseudo-code is shown in Figure 7. First, the machine learning model is generated from a precomputed initial dataset.Next, the objective function is predicted for a pseudo-input parameter set prepared separately using random numbers Latin hypercube sampling (LHS), and multiple input parameter sets that maximize the obtained values are extracted.The objective function is evaluated for the input parameter sets thus extracted and rigorously evaluated.The obtained calculation data are used again as training data to improve the accuracy of the machine learning model.In this way, the strategy is to increase the prediction accuracy in the region of the input parameter space where the objective function is large and updates the highest value.Since the evaluation of the edge of the data domain is important, we decided to use the Cubist algorithm [32], which is considered to be relatively resistant to extrapolation.

Benchmark problem 1: optimizations of analytic function
The algorithm implemented in MIOpt was verified on the minimization problem for the Rosenbrock function (Equation ( 1)).Here, the jobs were designed to perform four calculations in parallel at the same time.Seven cases shown in Table 1 were performed.The adopted j's were assumed to be logarithmically equally spaced {0.001, 0.005, 0.02, 0.1}.
The number of given data before starting optimization (initial data) was 10, and the number of BO trials was set to 30. Figure 8    points, which confirms the speed of optimization (the line with f ¼ 0:001 is shown in the figure for reference).In the sequential BO shown in Figure 8(a) and batch BO shown in Figure 8(b), the search was successful in the case of j(=0.001),whereas the search failed for other j values.The j-parallel BO (Figure 8(c)) was able to search widely around the optimal solution from Figure 8(c) (1) and to find the minimum value in an early iteration as shown in Figure 8(c) (2).In Cubist model update (Figure 8(d)), it can be seen that the search around the optimal value parameter was performed in Figure 8(d)(1); however, the point where the minimum value is shown in Figure 8(d)(2) was not reached.As a result, it was confirmed that j-parallel BO is effective for efficient optimization.

Benchmark problem 2: optimizations of welding conditions
In general, fracture of structures often occurs at the heat-affected zone (HAZ) of welds, and it is necessary to increase the strength of welds to ensure the integrity of structures.Therefore, it is important to evaluate the creep rupture time by the finite element method (FEM) with appropriate material constants, a heat source model, and a highly accurate fracture model, and by the identification of HAZ for heat input conditions.Various workflows have been developed for welds in MInt [33][34][35][36].On the other hand, accurate computational simulations are expensive and then a sufficient number of data points cannot be obtained.
Many studies have been carried out in the past to optimize welding conditions.Such studies include the optimization of heat source geometry parameters for GTAW (gas tungsten arc welding) in AISI1020 carbon steel and AISI304 stainless steel [37], the optimization of process parameters for 1/2" AA6061-T6 butt FSW (friction stir welding) joining [38], the optimization calculations of process parameters (current, voltage, root gap, and gas flow rate) for mechanical properties of TIG (tungsten inert gas) welding of stainless steel AISI 304 as base metal [39], the optimization based on the Taguchi method for the optimum welding parameters to obtain the ultimate tensile strength of the Ti-6Al-4 V welded specimens [40], and the optimization of geometric and welding process parameters of the gas metal arc welding process, using the design of the experiment and response surface method [41].The BO of the heat input conditions of the weld [42] has been carried out with several acquisition functions.

Workflow design and execution for creep rupture time evaluation
We previously developed a workflow for the weld creep rupture time of steel utilizing MInt [15] as shown in Figure 9.In this workflow, material constants, geometry, and heat source models are first determined, and thermal conductivity is analysed using the finite element method (FEM).Then, the maximum temperature for each element is obtained from the temperature history, and the HAZ, BM (base metal), and WM (weld metal) regions are identified by comparing the Ac 3 and Ac 1 temperatures.Next, creep rupture time is analyzed by FEM using the material constants corresponding to each of these regions, given boundary conditions such as load and temperature.The rupture time was evaluated using a timeconsumption law based on the principal stress with a resolution of 10 h [14].The FEM solvers used in the analysis are code_aster [43] for the thermal conductivity analysis and modified FrontISTR [44] for the creep rupture time analysis.
In this problem, the butt of the weld is I-shaped, as shown in Figure 10(a), and a three-pass butt weld is used.The heat source model is Goldak's double ellipsoidal distributed heat source model [45] with different Gaussian distributions from the center of the heat source to the front and back, as shown in Figure 10(b).The length of the heat source is divided 1:2 in the welding direction.
In the FEM calculations, a symmetrical part (1/2 model) from the center of the weld was analysed, and 32 cores per job were used to perform parallel calculations.To evaluate the performance of the algorithms proposed in this study, the problem of optimizing the heat input conditions to maximize the creep rupture time at the weld is verified.

Optimization and the largest creep rupture time
Optimization calculations were performed to find the heat input conditions that maximize the creep rupture time.The input parameters are the heat input conditions during welding (heat input, and width and depth of each layer).In this study, the bead width and thickness were fixed at 4 mm and 8 mm, respectively.Therefore, the heat input depth of the third layer was determined by the heat input depths of the first and second layers; therefore, it was not treated as an independent variable.Therefore, 8 variables were used, and search ranges are shown in Table 2.
Let x i be the input variables (heat conditions) and χ be the input parameter space.The problem of maximizing (T r ) the creep rupture time (T max r ) can be written as, T max r :¼ arg max The applied stress was 100 MPa and the test temperature was 873 K in the calculations.
In the thermal conductivity analysis, the material constants for the entire model (elastic constants (E 0 ), Poisson's ratio (ν), phase transformation temperatures (Ac1, Ac3), the coefficient of thermal expansion (CTE), and thermal conductivity (λ)) shown in Table 3 were used.
The strain rate during creep can be written as a function of stress σ and time t shown as Equation ( 4) using the Norton-Arrhenius type formula.where T is temperature, R is the gas constant, Q � is the apparent activation energy for the minimum creep rate, n � is the stress index, and the parameter A � is the material constant.The material for the analysis is assumed to be 2.25Cr-1.0Mosteel as shown in Table 4.

_ εðσ; TÞ
As described in Section 4.1, we used four algorithms, sequential BO, batch BO, j-parallel BO, and Cubist model update, for the optimization calculations and performed seven different calculations shown in Table 1.The initial data set was calculated using 49 data generated from an 8-variable, 7-level L49 orthogonal table (orthogonal array, OA) on the basis of experimental design methods to ensure uniform data sampling for the input parameter space.The largest value of rupture time in the preliminary L49 evaluations was 580 h.The number of iterations for optimization was set at 12, assuming that the actual time required for optimization is approximately two weeks.

Comparison of results by algorithm
The updated histories of the highest values for various optimization calculations are shown in Figure 11 [(a) sequential BO, (b) batch BO(j = 0.001), (c) batch BO(j = 0.005), (d) batch BO(j = 0.02), (e) batch BO(j = 0.1), (f) j parallel BO, and (g) Cubist model update].'The number of iteration' on the horizontal axis shows jobs (four in this case) running simultaneously in the same iteration.Therefore, the iteration can be taken as the actual time required.
For sequential BO (Figure 11(a)) and batch BO (Figure 11(b)-(e)), using j = 0.02, we looked for input parameter sets that achieve rupture times exceeding 600 [h].For j = 0.001 or j = 0.005, we were unable to find a condition that achieves a rupture time exceeding 600 [h].This means that knowing the appropriate j is very important for computational efficiency.On the other hand, the j-parallel BO (Figure 11(f)) obtained a rupture time of 640 [h] as the highest value at iteration number 7, which is 10% higher than the highest value (580 [h]) obtained under the L49 condition.We found condition with a rupture time of 640 [h] by the Cubist model update (Figure 11(g)) algorithm.The difference between a calculation that performs multiple sequential simultaneously and a j-parallel BO is whether or not the data are synthesized in each iteration.Comparing sequential BO (Figure 11(a)) and j-parallel BO (Figure 11(f)) shows that the search is more efficient through the procedure of data synthesis at each iteration.This j-parallel BO was evaluated several times while changing the initial values, and a similar convergence situation was obtained.
Next, principal component analysis (PCA) was performed on all input data used to calculate rupture time.Among them, the results of the PCA of the initial data (L49) and each j used in the sequential BO are shown in Figure 12 [(a) j = 0.001, (b) j = 0.005, (c) j = 0.02, and (d) j = 0.1].The PC1 and PC2 axes have an explanatory performance of 36%, and the dispersion of the input values is considered on a plane composed of these two axes.Since the initial data L49 are selected uniformly in the design of the experiment, it has a large dispersion with respect to the first and second axes of the principal components.In sequential BO, the search range varies depending on the value of each j, but it is not wide for the spread of L49.However, depending on j, it is confirmed to have a wide range in the input space, and it can be considered that a wide range of input space can be searched.Therefore, the choice of j has a significant impact on the efficiency of optimization.
Figure 13 shows the results of PCA of data for batch BO [(a) j = 0.001, (b) j = 0.005, (c) j = 0.02, and (d) j = 0.1].The results depend on the value of j, the same as those of sequential BO.Therefore, to improve the efficiency of the optimization, the optimal j must be prepared in advance.
Finally, we review the parameters of the j-parallel BO and Cubist model update proposed in this study.Figure 14 shows (a) j-parallel BO, (b) Cubist model update.It can be judged that j-parallel BO can search  a wide range of input space.The number of data points can be increased by synthesizing data at each iteration, and as shown in Figure 12, it is possible to find parameter spaces that cannot be searched by utilizing a single j.Thus, when performing optimization in calculations where it takes time to acquire a single data point, we found that the strategy of increasing the number of data points by maximizing the use of computational resources while maintaining a balance between exploration and exploitation works well.On the other hand, when using the Cubist model update shown in Figure 14(b), it can be seen that the search is around the parameter space that achieves the highest value.When using the Cubist model update, the highest value could not be obtained in the end.
As a result, it was confirmed that performing BO while synthesizing data, such as j-parallel BO, is more effective than multiple individual sequential BOs or batch BO with fixed j (Appendix A).

Data
The calculated data are stored on GitHub, including data for the Rosenbrock function optimization and the weld heat input condition optimization data containing L49 calculation.https://github.com/materialsintegration/MIOpt/tree/main/opt_heatcon

Summary
Accurately solving forward problems in structural materials requires the manipulation of a large number of physical models and many input -output descriptors.In the forward problem analysis infrastructure (MInt), these are described and implemented as a workflow.In this study, we have developed an environment (MIOpt) that allows them to link as optimization problems.Therefore, once the forward problem is designed, the backward problem can be carried out immediately.MInt/MIOpt can be used free of charge for researchers at universities or public research institutes.
Furthermore, we have implemented our own optimization algorithms, j-parallel BO and Cubist model update.In particular, j-parallel BO was improved to maintain a balance between exploration and exploitation and was found to improve optimization efficiency.The problem used as an example in this study is the minimization of an analytical function and the verification of the problem of optimizing the heat input conditions to maximize the weld creep rupture time.The results show that j-parallel BO is more efficient in finding the optimal set of input parameters than to running sequential BO simultaneously.
Although it is not possible to generally determine the value of j that works for every problem, we believe that the method of using multiple j's and synthesizing data works effectively.

Figure A2
. Geometrical features of weld shape that achieve more than 610 h of rupture time.In several algorithms, the parameter set was not found; in such a case it is indicated as "empty".

Figure 1 .
Figure 1.Schematic of workflow for forward problem in "MInt".The blue boxes indicate the input ports (x 1 , x 2 ), the gray box the output port (f ), and the yellow box the prediction model (function).

Figure 2 .
Figure 2. Linking MIOpt and MInt through MInt-API to perform backward problem.

Figure 3 .
Figure 3. Image of the linkage between MIn and MIOpt.A forward problem workflow set in MInt is connected to MIOpt to solve the backward problem.In MIOpt, the workflow selection, initial data, and constraint conditions that are the setting information of the backward problem are possible.

Figure 5 .
Figure 5. Schematic of selecting next candidate point with an acquisition function (expected improvement, EI) based on Gaussian process.
shows the results of the search using each algorithm [(a) sequential BO, (b) batch BO, (c) j-parallel BO, and (d) Cubist model update].'(1)' in the left panels of Figure 8(a-c) are sampled points, and '(2)' in the right panels of the figures indicates the evolution of the evaluated values for the sampled

Figure 7 .
Figure 7. Algorithm for Cubist model update optimization.

Figure 9 . 10 .
Figure 9. Workflow for creep rupture time analysis implemented in MInt.

Figure 12 .
Figure 12. Results of extracting initial data L49 and sequential BO data for each j from PCA of all data [(a) j = 0.001, (b) j = 0.005, (c) j = 0.02, and (d) j = 0.1].Black circles are the L49 data, red circles the highest points in each sequential BO, red open circles the highest in each BO.

Figure 13 .
Figure 13.PCA, of L49 and several batch BOs [(a) j = 0.001, (b) j = 0.005, (c) j = 0.02, and (d) j = 0.1].Black circles are the L49 data, orange circles the highest points in each batch BO, red open circles the highest points in each BO.

Figure 14 .
Figure 14.PCA of j-parallel BO and Cubist model update [(a) j-parallel BO, and (b) Cubist model update].Black circles are the L49 data, green or blue circles the highest points, and red open circles the highest in each optimization.

Table 1 .
Calculation cases to verify the algorithms.

Table 2 .
Input variables, survey ranges, and descriptions.

Table 3 .
Material properties of 2.25Cr-1Mo steel for thermal conductivity analysis.