Crack Identification by Digital Image Correlation Method Using Crack Shape as Prior Information

The integrity of large electrical machines is regularly inspected for cracks to verify their integrity [1]. Cracks in electrical machines occur in both visible and invisible locations. One of the electrical machines is a turbine generator. Cracks can occur in the retaining ring of turbine generators [2]. The retaining ring holds the centrifugal deformation of the rotor coil. The retaining rings are cylindrical, with a thickness of about 10 to 20-mm and a diameter of about 1-m. Cracks in the retaining ring occur where it is shrink-fitted to the rotor [2] and are not easy to visually inspect. To detect cracks in invisible locations, an electrical machine is shut down for a few weeks, disassembled, and inspected [3]. If cracks in invisible locations could be inspected without disassembling electrical machines, shutdown periods would be reduced.

In recent years, robotic inspections have increased to avoid such disassembly. Robots enter narrow spaces inaccessible to humans [3–6], and the inspection devices mounted on them are as thin and small as possible. For example, the hammering inspection device for robots entering turbine generators [7] is a thickness of approximately 20 mm or less [8].

Nondestructive inspection methods for cracks in invisible locations include hammering inspection [9], ultrasonic inspection [10], infrared inspection [11], microwave inspection [12], X-ray inspection [13], eddy current testing [14], and piezo-electric sensor [15]. Hammering and ultrasonic inspections use sound waves; infrared, microwave, and X-ray inspections use electromagnetic waves. Piezo-electric sensor vibrates the inspected object by applying alternating current. Eddy current testing uses the changes in an electromagnetic field to identify cracks. All of these inspection methods are changed by waves or electromagnetic fields from cracks. A smaller generating device results in weaker waves or electromagnetic field input. A smaller receiving device is less accurate for measuring waves or electromagnetic fields. Achieving downsizing and sufficient accuracy is difficult for both the generating and receiving devices.

In these nondestructive inspection devices, a device maintains contact with an object to reduce input losses. An inspection device might get caught in a narrow space. If the narrow space is several meters deep, removing the robot is complicated.

Camera-based visual inspection is a noncontact inspection method. A camera can be easily miniaturized for finding visible cracks [16,17], which must be predicted based on their image changes. Digital holographic microscopy(DHM) uses light interference to precisely measure surface deformation [18]. The DHM measurement device needs an optical interference system. The DHM measurement device is not easy to miniaturize. A scheme that quantifies image changes is the digital image correlation (DIC) method [19], which identifies deformations from images. It determines the amount of deformation based on the changes in a pattern of a photographed surface. Several methods for predicting cracks in invisible locations use the DIC method. One of the conventional methods is to prepare deformations due to cracks with numerical analysis in advance [20]. This method compared measured deformations with the deformation of each crack's shape. The estimated result is the one closest to the measured deformation from the deformation for each crack's shape. With this method, it is difficult to estimate cracks that are not present in advance data or when the measured deformations have noises.

The other conventional method is to solve the inverse problem. The inverse problem estimates the displacement of crack surfaces from surface deformation. Since measured deformations contain noise, the inverse problem has an ill-condition. A method for solving the inverse problem with an ill-condition is to predict crack shapes using topology optimization [21–23]. Such topology optimization research leverages full-field response data obtained by DIC in a topology optimization framework to reconstruct the internal damage in structural members [21]. The authors extended their research to demonstrate the feasibility and performance of their proposed method for a set of large-scale structural steel beams with/without buried defects using a full-field 3D DIC sensor approach [23]. This conventional method suffers from the high computational costs of repeatedly computing a finite element method in the optimization.

To overcome this ill-condition, Amaya et al. [24,25] proposed regularization that considers the physical relationship between the surface deformation and the crack deformation and the physical constraint between the displacement and the force at a crack and its perimeter (ligament). This regularization method is based on the Joint Estimation Maximum a Posteriori (JE-MAP) method [26,27], which was developed to improve the identification accuracy of body tissues by X-ray CT results. The JE-MAP algorithm alternates between Maximum a Posteriori (MAP) estimation [28], which has physical constraints as the prior information, and the grab-cut (GC) method [29], to avoid ill-conditions. In this paper, the regularization method that developed the JE-MAP method is called the JE-MAP method. However, depending on the distribution of noise contained in the measured deformations, cracks were sometimes estimated at locations inconsistent with fracture mechanics.

Therefore, in this study, we developed a method that identifies cracks in invisible locations by combining Expectation a Posteriori (EAP) estimation [30,31] and a DIC method. We consider the inverse problem that identifies crack shape from observation of the surface deformations by DIC measurements. In order to regularize the inverse problem, EAP with prior regarding crack shape candidates is proposed. The crack shape candidates are generated with numerical analysis so that these candidate crack shapes may follow the fracture mechanics. This regularization method is called crack shapes prior EAP(CSP-EAP). Since CSP-EAP is a random variable for both the candidate shapes and the measured deformations, this method can identify even cracks that are not present in crack shape candidates or when the measured deformations have errors. In order to demonstrate the effectiveness of the proposed method, the validation test, which compares with the conventional L1-norm regularization method and the JE-MAP method, is performed. The specimens and test conditions in this study are based on cracks in the retaining rings of turbine generators and the strain due to shrink-fit [2,32].

The inverse problem in this study is the same as the problem solved by our previous study [25]. Figure 1 shows an overview of the inverse problem.

The inverse problem estimates the unknown data from the observed data in Fig. 1. The observed data $εZ$ are the Z-direction strain changes of the observed surface. The unknown data are Z-direction displacement u of the virtual plane. The observed data $εZ$ is the change with the presence/absence of a crack under tensile loading (Fig. 1).

Here, $β$ is the error term involving the numerical analysis error and the measurement error. $εZ˜, u$⁠, and $β$ in bold are discretized vectors. $H$⁠, obtained by the finite element method or other methods, is a constant matrix representing the relationship between $εZ$ and $u$ without $β$⁠.

In this study, our proposed regularization method uses the EAP method with prior information on the crack shapes. The prior information is the candidate crack shapes and the observed surface deformations from them. They are generated by numerical analysis. The candidate crack shapes must be consistent with fracture mechanics. The posterior distribution is obtained from the observed data and the prior information. The estimated displacement of the virtual plane is the expected value of the posterior distribution. Here, we describe the computational flow for our proposed regularization method, the creation of candidate geometries, and the determination of hyperparameters for the regularization method.

In this section, we describe a numerical analysis method for creating candidate crack shapes. Several candidate shapes can be created consistently with fracture mechanics. When we perform the crack propagation simulation, various fracture criteria can be employed, such as: the stress based criterion [33], the COD criterion [33], the criterion with stress intensity factor [34], with J-integral [35], with plastic strain [36], or inelastic strain energy density [37,38]. Also, several crack propagation techniques are available, such as X-FEM [39], delete and recreating mesh approach [40–42]. Based on these numerical simulation techniques, usually the information of material properties [43] are taken into account.

In this study, the force on the virtual plane is used as a parameter to determine the crack propagation. The crack's region is determined by a set of conditions based on this force. Candidate shapes are obtained by modeling the object using the finite element method. The virtual plane's crack shape is reflected as a boundary condition in the finite element method. Force vector f, displacement vector u, and Z-directional strain variation $εZ$ of the observed surface are obtained. Perturbation analysis [33,44] and other methods create candidate geometry to reduce the computational cost. Latent variables, which are introduced to indicate the candidate crack shapes, are indicated by 0 and 1 for the crack and ligament regions on the virtual plane. The boundary conditions of the virtual plane are determined based on the latent variables. Figure 2 shows the flow of creating candidate shapes. The starting point of the candidate shapes is set in the virtual plane. The latent variable at the starting point is 1, and 0 otherwise. The boundary conditions of the virtual plane are determined by the latent variables. The displacement, the virtual plane's load, and the strain distribution of the observed surface are obtained by numerical analysis with these boundary conditions. Based on the force distribution on the virtual plane and predefined conditions, the next range of cracks is determined. The latent variables are set to 1 for that range. This process is repeated to create candidate shapes.

Prior information includes the observed surface strains, the latent variables, the displacements, and the forces of the virtual plane. The EAP method estimates the expected value of a posterior distribution, which is defined from the probability distributions of the observed surface strain and the surface strains for the candidate shapes. The EAP method is identical to the MMSE method [45,46] when the likelihood distribution is a normal distribution with constant variance and zero covariance, and the prior distribution is vague.

In this study, a parameter is number v that indicates the candidate crack shapes. First, the likelihood distribution is obtained of observed data $εZ˜$ and parameter v. The likelihood distribution is a probability distribution of the difference between observed strain vector $ε(v)$ of the candidate shape indicated by parameter v and measured strain vector $εZ˜$⁠. Figure 3 shows an overview of the likelihood distribution of the prior information and the observed data.

$Catv$ is a probability density function of a categorical distribution. Since the probability is identical for all v, the prior distribution is a discrete uniform distribution.

Another estimated method is where v is selected for the greatest similarity with the observed data. In the proposed regularization method, the expected value of the candidate shapes for each discretized point can be predicted by combining multiple candidate shapes, even for crack shapes that do not exist in the candidate shapes.

Figure 4 shows the flow of the CSP-EAP. Latent variable z determines the starting point of the candidate crack shape. The displacement, the force of the virtual plane, and the observed surface strain are obtained under the boundary conditions determined by the latent variables. The crack propagation is repeated based on the obtained forces. When a crack penetrates the virtual plane, the creation of candidate crack shapes is completed. The product is calculated of the likelihood and prior distributions. A marginal distribution of the observed data is summed over the product of the likelihood and prior distributions. The posterior probability is calculated from the product of three distributions: likelihood, prior, and marginal. The expected value of displacement $u˜$ is calculated from the posterior probability.

where Num is the number of discretized data points on the observed surface.

$σ2$ is a value that minimizes MSE. The method for determining $σ2$ when calculating the likelihood distribution is shown in Sec. 3.3.

In this study, the specimens and the loading conditions are the same as the problem solved by our previous study [25]. Figure 5 shows an illustrative diagram of a 50-mm wide, 24-mm thick flat plate. The identified crack exists on the virtual plane, indicated by the dashed line at the plate's center (Fig. 5). In our study, we measured the strain changes of the observed surface in Fig. 5 without cracks with the DIC method. A load of 247-kN is applied to the specimen in Fig. 5 to generate a strain of 1000 μ in the Z-direction.

The following section describes a numerical model for creating the candidate crack shapes shown in Sec. 2.2.1 and the conditions for determining the extent of the crack propagation from the force distribution on the virtual plane.

Figure 6 is the finite element model from which the candidate crack shapes were created. The finite element model is the same as the problem solved by our previous study [25]. The vectors of strain $εZ$⁠, displacement $u$⁠, and force $f$ in Fig. 6 are the values at the nodes of the finite element method. The virtual plane in Fig. 6, the width of 50-mm and the height of 24-mm, is divided into $n=$25 and $m=$12 subdivisions. The observed surface in Fig. 6, the width of 50-mm and the height of 24-mm, is divided into $n=$25 and $p=$15 subdivisions. The flat plate material is quenched and tempered of S45C steel (C45 steel: ISO) with a proof stress of 490 MPa or higher. Young's modulus and Poisson's ratio of the plate are set to 206 GPa and 0.3. The finite element method solver is ansys 19.2 [47].

Figures 7 and 8 are schematic diagrams of the conditions of this paper. Figure 7 shows two cases: a point with the greatest nodal force and another with the points with the greatest and second greatest nodal force. The latent variable determines the crack initiation point. The displacement, the force of the virtual plane, and the observed surface strain distribution are determined under the boundary conditions containing the crack initiation point. The force distribution is converted to the nodal force change rate using Eq. (14). The two points with the greatest nodal force change rate are the next cracks, and the latent variable is set to 1, which is the crack. Figure 8 shows the case where the nodal force change rate is greater than 90%. Here, the crack initiation points are two consecutive points in the X-direction of the virtual plane. As in Fig. 7, the force change rate is determined, and the conditions for propagating the crack were the point with all nodes with a nodal force change of 90% or greater. The Same is true for a nodal force change of 99% or greater, which is a smaller range from which to choose. This process is repeated until the crack penetrates in the Y-direction of the virtual plane.

The candidate crack shapes for Sec. 2.2.1 are the result of eliminating the same crack shape among the candidate shapes obtained by each condition. Unlike the JE-MAP method, our method does not estimate a crack shape at a location inconsistent with fracture mechanics, since it is obtained by adding up the candidate crack shapes that are consistent with fracture mechanics. An example of a candidate shape used in this paper is shown in Fig. 9. Candidate shapes are shown in white for cracks and in black for noncracks. The total number of candidate shapes is 4028. The crack propagation based on a J-integral or a stress intensity factor may be used to generate candidate shapes.

In this section, the parameters in Sec. 2.2.3 are determined. The parameter is a value that minimizes the MSE in Eq. (13). A range is predetermined within which the parameters are searched. Standard deviation σ corresponds to the error of the strain distribution. σ is specified by considering both the numerical analysis error and the measurement error. The numerical analysis error is assumed to be 1% at most [48]. We assume σ in the range of 1 μ to 100 μ. The parameter is optimized using Bayesian optimization [49].

The shape of the specimens is the same as the problem solved by our previous study [25]. Figure 10 shows the shape of the specimens which were measured by the DIC method. A specimen without a crack is shown in Fig. 10(a), and a specimen with one is shown in Fig. 10(b). The specimens include the part gripped by the testing machine. The crack specimen has a semi-elliptical slit that is 0.25-mm thick, 24-mm wide, and 9.8-mm depth. A semi-elliptical, longitudinal slit was made in its center by electrical-discharge machining. Although the stress state at the crack tip differs between the slit and the crack, the strain distribution on the observed surface is almost the same. Because the observed surface is away from the slit tip. The situation measured by the DIC method with a tensile load of 247 kN is shown in Appendix  A.1. This loading condition causes strain in the specimen equivalent to the strain due to shrink-fit in the retaining ring of the turbine generator [2,32] assumed in the verification of this paper. The specimen surface pattern for measurement by the DIC method, the imaging equipment, and the software settings for the DIC method are the same as those used in our previous study [25].

The observed data for the inverse problem are the strain changes obtained by the DIC method from the images before and after the crack initiation, similar to the validation of the JE-MAP method [25]. The observed data is the amount of strain change in the Z-direction with and without the crack. Figure 11 shows the difference in strain in the Z-direction between the case without a crack and the case with a crack, which is the observed data. The strain distribution in Fig. 11 is 0 in the Z-direction at the lower left and 0 in the X-direction at the left end of the specimen width. Unlike the JE-MAP method [25], the observed data is within the strain distribution shown by the dashed lines in Fig. 11. The observed data are the discrete values of the strain distribution shown in Fig. 11 divided into the element sizes shown in Fig. 6.

Figure 12 shows the strain distribution used as the observed data in the inverse analysis. The method for determining the observed data is the same as in the JE-MAP method [25] and is shown in Appendix  A2. Since the strain distribution is measured by the DIC method, the strain distribution at the specimen edge in the X-direction is missing. The Z-position of the observed data is defined as the distance from the 32 mm Z-coordinate in Fig. 11. The observed data in the X-direction are the area without the missing area measured by the DIC method. $H$ was recalculated based on the range of observed data.

In this section, we solve the observed data shown in Fig. 12 using a generic method: L1-norm regularization. Details of the L1 norm regularization are given in Appendix  B. The regularization parameter λ was set to $1.00×10−6$⁠. Figure 13 shows the estimated results of the z-direction displacement of the virtual plane. The 0 of the Y-coordinate in Fig. 13 is the surface where the crack opens, and the 0 of the X-coordinate in Fig. 13 is the left edge of the virtual plane of the 1/2 symmetry model in Fig. 6. The Z displacement of the virtual plane is positive for the cracked part and zero for the uncracked part. Figure 14 shows the crack extracted from Fig. 13. The extracted crack is the area of displacement greater than 10% of the maximum displacement of the virtual plane and is indicated by 1. The area other than the extracted crack is indicated by 0. The zero points of x and y in Fig. 14 are the same as in Fig. 13. Figure 15 shows the ground truth displacement distribution obtained using the finite element method on a model with the same crack as the test specimen. The zero points of x and y in Fig. 15 are the same as those in Fig. 13. The L1-norm regularization estimates two cracks smaller than in Fig. 15. Unlike Fig. 15, the extracted cracks are away from the free edge of the specimen. As in our previous study [25], the L1-norm regularization did not take into account the error included in the strain changes of the observed data, causing a result different from the ground truth.

In this section, we solve the observed data using the JE-MAP method [25]. The hyperparameters of the JE-MAP method are zConst, $σw2$⁠, and the initial value of μ_uc. Table 1 shows the values of the hyperparameters. The hyperparameters in Table 1 were obtained by Bayesian optimization as in our previous study [25]. Figure 16 shows the estimated results of the z-direction displacement of the virtual plane. The x- and y-axes and the direction of displacement indicating cracks in Fig. 16 are the same as in Fig. 13. Figure 17 shows the distribution of the latent variable in the prior distribution used to compute Fig. 16. The latent variables are 1 for the part indicating the crack and 0 for the rest of the crack, as in Sec. 2.2.1. The x and y zero points in Fig. 17 are the same as in Fig. 13.

From Fig. 16, positive displacements are predicted in areas that are predicted cracks, and zero or negative displacements are predicted in areas that are not. Compared to our previous study [25], there is a larger error in the distribution of the estimated displacements. Even the cracks shown by the latent variables in Fig. 17 predict multiple cracks, which is very different from the cracks in Fig. 15.

Figure 18 shows the distribution of the rate of the force change on the hypothetical surface for the crack shape in Fig. 17. There are areas around the crack where the rate of force change is small, predicting a crack shape that is not following fracture mechanics. In the previous sturdy [25], cracks were also estimated cracks from DIC measurements, but it was found that the accuracy of crack estimation was reduced when the measurement errors were different. The CSP-EAP method shown in Sec. 2.2 is a regularization method that can estimate cracks according to fracture mechanics, taking into account measurement errors in the observed data.

In this section, we solve the observed data using the CSP-EAP. Table 2 shows the value of σ for the parameter presented in Sec. 3.3. The Bayesian optimization that obtains the value of σ is shown in Appendix  C. σ is specified by considering both the numerical analysis error and the measurement error. The standard deviation of the error σ is close to the standard deviation of the strain distribution measured by DIC, $2.6×10−5$ (Appendix  A.2).

Figure 19 shows the estimated results of the z-direction displacement of the virtual plane. Figure 20 shows the crack extracted from Fig. 19. The crack area is defined as Fig. 14 shown in Sec. 5.1. The maximum point of the displacement distribution in Fig. 19 is close to the ground truth in Fig. 15. As shown in Fig. 20, only one crack is estimated, and most of the area estimated to be a crack is included in the ground truth in Fig. 15.

The identified size of the crack is compared to the ground truth. The area was determined from the number of elements in the cracked regions shown in Figs. 20 and 15. One element in the distribution figures in Figs. 15–20 has a width of 2 mm and a height of 2 mm. The area identified by the CSP-EAP is 148 mm². The ground truth shown in Fig. 15 is 224 mm². The area identified by the CSP-EAP is 66% of the ground truth. The sizes identified by the CSP-EAP are a maximum crack width of 24 mm and a maximum crack depth of 10 mm. The maximum width and depth were identified to be 90% and 83% of the ground truth in Fig. 15. The ratio of the estimated crack the maximum depth to the maximum width is in 92% agreement with the ground truth. Therefore, the CSP-EAP can identify a crack's size and location from the observed data calculated by the DIC method.

Figure 21 is the distribution of the rate of the force change on the hypothetical surface for the crack geometry in Fig. 20. The zero points of x and y in Fig. 21 are identical as in Fig. 13. As shown in Fig. 21, the force change around the crack is large. As shown in Fig. 18, the JE-MAP method sometimes identified cracks even at locations where the force change around the crack was not large. The CSP-EAP does not cause the problems that occurred in the JE-MAP method.

We believe that the following factors contributed to the improved estimation of the CSP-EAP. The CSP-EAP is the expected value of the posterior distribution of the candidate shapes under the conditions of the observed data. Many candidate shapes were prepared by speeding up their creation by propagating cracks based on the force at the ligament. Among the many candidate shapes, we included a candidate near the crack that produces the observed data. The error in the observed data was introduced as a variance in obtaining the posterior distribution of the candidate shapes under the conditions of the observed data. The posterior distribution can be obtained by taking into account the errors in the observed data. Its expected value for each point of the discretized virtual plane can estimate a crack's shape that is not in the candidate shapes.

In this study, we developed a method that identifies cracks in invisible locations by combining EAP estimation and a DIC method. The following are our results:

• Our proposed method identified the shape of cracks in invisible locations by obtaining prior information of crack shapes prepared a priori by numerical analysis and the surface deformation obtained by the DIC method.

• The proposed scheme is the EAP method, which calculates the expected value of the posterior distribution of the prior information of a crack's shape relative to the measured surface deformation.

• The candidate crack shape was generated by propagating cracks based on the load distribution of the ligament, and prior information was generated with low computational cost.

• The effectiveness of the proposed method was compared and validated using the specimen simulating turbine generator components and loading conditions.

• The CSP-EAP estimated the displacement on the crack-containing plane more accurately than L1-norm regularization and the JE-MAP method in the inverse problem where the observed data were strain distributions measured by the DIC method with the measurement error.

• The CSP-EAP estimates cracks more accurately than the L1-norm regularization method and the JE-MAP method for the inverse problem where the observed data were the strain distribution measured by the DIC method.

• The CSP-EAP identified the crack size and location based on surface deformation, including the measurement error, which the JE-MAP method failed to identify.

The JE-MAP method can be used in combination with the CSP-EAP. For example, the JE-MAP method can be used to estimate the number and the approximate shape of cracks. Based on the JE-MAP results, candidate shapes are generated and estimated by the CSP-EAP. By narrowing down the candidate shapes to be used for CSP-EAP, it may be possible to estimate the crack shape even if the measurement error is large. In the future, the validity of the combination method of the CSP-EAP with the JE-MAP will be verified for cracks of various shapes, loading conditions, and specimen geometry, e.g., multiple cracks, and ring specimen.

DIC =

digital image correlation

JE-MAP =

joint estimation maximum a Posteriori

EAP =

expectation a Posteriori

CSP-EAP =

crack shapes prior EAP

$L$ =

likelihood distribution

$Pri$ =

prior distribution

$Pos$ =

posterior distribution

$εZ$ =

strain changes in the Z-direction of the observed surface

u =

displacement of the virtual plane in the Z-direction

$εZ˜$ =

vector of the measured strain changes in the Z-direction of the observed surface

$u$ =

displacement vector of the virtual plane in the Z-direction

$β$ =

discretized vector of the measurement errors

$H$ =

constant matrix representing relationship between $εZ$ and $u$

v =

number indicating the candidate shape

$v$ =

vector of numbers indicating candidate shapes

$u˜$ =

crack surface displacement maximizing the posterior probability shown in Eq. (8)

$σ2$ =

variance of measurement error

$f$ =

discretized Z-direction force vector of the virtual plane

MSE =

mean square error

λ =

regularization parameter of L1-norm regularization

Figure 22 shows the situation caused by the tensile 247-kN load where the strain is measured by the DIC method. Figure 23 shows the strain distribution in the Z-direction of the observed surface calculated by the DIC method at a load of 247 kN. Figure 23(a) shows the specimen's strain distribution without a crack, and Fig. 23(b) shows it with a crack. The no-crack specimen has a strain distribution without any characteristic distribution. The strain distribution had a mean value of 1,005 μ and a standard deviation of 26 μ. The specimens with cracks showed a decrease in strain near the strain distribution's center. This decrease was caused by a crack on the opposite side of the observed surface.

Figure 24 shows the method that determined the crack's location in the Z-direction based on the strain changes in Fig. 23. Figure 24(a) shows the method that determines the Z-coordinate where the strain changes are maximum. The discretized strain change was divided into each X-coordinate. The maximum Z-coordinate is obtained for each divided strain change. Figure 24(b) shows the frequency distribution of the Z-coordinates with a maximum strain change. The virtual plane is located at the Z-coordinate with the highest frequency of maximum strain change. The strain change's distribution is divided by the virtual plane into positive and negative directions of the Z-axis.

λ denotes a regularization parameter. Regularization parameter λ was calculated by k-fold cross-validation with k = 5. Appendix Eq. (15) was solved using the Fast Iterative Shrinkage-Thresholding Algorithm (FISTA) method [51]. The method for determining the regularization parameter λ is shown. The observed data were partitioned by cv partition, a data creation function for cross-validation [52] in matlab.

λ was searched for in the range of $1×10−8$ to $1×10−5$⁠. One hundred λ were extracted from the search range at logarithmically evenly spaces. Cross-validation was performed for all the extracted λ. Appendix Fig. 25 shows the mean and standard deviation of MSE calculated by cross-validation. The vertical axis shows MSE and the horizontal axis shows λ. In Appendix Fig. 25, the mean and standard deviation of the MSE calculated by cross-validation for each value of λ are shown as circles and error bars. λ was determined to be $1.00×10−6$ by applying” The One Standard Error Rule” to the results in Appendix Fig. 25.

The σ of the CSP-EAP and the parameters of JE-MAP were determined by the Bayesian optimization [53] of matlab function bayesopt. The following are the bayesopt arguments that were changed from the default values:

• Acquisition function: Uses expected-improvement and enables modification of acquisition function to escape a local objective function minimum.

  ”AcquisitionFunctionName” =“expected-improvement-plus”

• Specify deterministic objective function: true

  (Objective function is specify-deterministic.)

• Objective function evaluation's limit: 60