From specified product tolerance to acceptable material and process scatter: an inverse robust optimization approach

Production efficiency in metal forming processes can be improved by implementing robust optimization. In a robust optimization method, the material and process scatter are taken into account to predict and to minimize the product variability around the target mean. For this purpose, the scatter of input parameters are propagated to predict the product variability. Consequently, a design setting is selected at which product variation due to input scatter is minimized. If the minimum product variation is still higher than the specific tolerance, then the input noise must be adjusted accordingly. For example this means that materials with a tighter specification must be ordered, which often results in additional costs. In this article, an inverse robust optimization approach is presented to tailor the variation of material and process noise parameters based on the specified product tolerance. Both robust optimization and tailoring of material and process scatter are performed on the metamodel of an automotive part. Although the robust optimization method facilitates finding a design setting at which the product to product variation is minimized, the tighter product tolerance is only achievable by requiring less scatter of noise parameters. It is shown that the presented inverse approach is able to predict the required adjustment for each noise parameter to obtain the specified product tolerance. Additionally, the developed method can equally be used to relax material specifications and thus obtain the same product tolerance, ultimately resulting in a cheaper process. A strategy for updating the metamodel on a wider (noise) base is presented and implemented to obtain a larger noise scatter while maintaining the same product tolerance.


Introduction
Every metal forming process is influenced by the scatter of both material and process conditions. Robust optimization is an essential approach to enhance the production quality in forming processes in the presence of the scatter. Recently, robust optimization has been implemented for cold roll forming [1], stretch-drawing of a hemispherical cup [2], extrusion-forging [3], stamping process [4], and V-bending O. Nejadseyfi o.nejadseyfi@utwente.nl [5]. The goal in robust optimization is to find a design setting which reduces the variability of the product in the presence of noise, while a specified percentage of production satisfies constraints.
In this approach, the input parameters to a process simulation are separated in two categories, design parameters (x) and noise parameters (z). Design parameters can be adjusted during the forming process while noise parameters are difficult to control or are out of control. The noise parameters can be either from the material or from the process. For example, the thickness variation or the scatter in hardening behavior of the material are typical material noise whereas the variation in forming temperature or lubrication can be considered as process noise. The measured or estimated input noise is then propagated via the process and variation of the output is calculated. For this purpose, methods such as Monte Carlo (MC) and its variations, Taylor-series expansion, Gaussian Quadrature (GQ), polynomial chaos, stochastic collocation, and analytical methods have been widely used [6]. An optimization algorithm can be employed to find the optimum design at which a minimal product variation in the presence of input noise is obtained. A common approach is to minimize the variation while setting the mean to target simultaneously [5].
A tighter product tolerance is only achievable by requiring less scatter of noise parameters. Finding a solution to reduce the scatter of input noise, and implementing those solutions usually incurs additional costs [7]. For example, for a smaller variation of thickness, it is necessary to precisely monitor and control sheet production, and coils that do not satisfy the requirements will be excluded from delivery or will be sold to less stringent clients at a lower price.
The subject of this research is to determine the acceptable material and process scatter from a specified product tolerance by inverting the robust optimization method. This inverse problem is regarded as the reverse process of a robust optimization problem and is referred to as tailoring of scatter. The non-linearity of the process, the presence of multiple noise variables and the interactions between them, are the main challenges in inverse approach that must be addressed.
To solve inverse problems in forming-related processes, various methods have been implemented; for example, straightforward but computationally expensive brute force strategy [8], Genetic Algorithm (GA) for global optimization [9], and gradient based optimization methods [10]. Most of those studies consider deterministic inverse approaches; meaning that the effects of variation of input and output are ignored.
In this article a gradient-based optimization algorithm, sequential quadratic programing (SQP), is implemented to solve the inverse robust optimization problem in Matlab. The gradients of objective function and constraint are evaluated analytically as described in [11]. An inverse approach is proposed and is applied to forming a lab-scale automotive part. A metamodel of this process is made for rapid calculation of uncertainty propagation as is common practice when dealing with expensive computer simulations. For example, Keating et al. [12] implemented the Null-Space Monte Carlo (NSMC) method on a meta-model of a large non-linear problem for parameter estimation and uncertainty propagation. Pacheco et al. [13] used MC integration on the metamodel built from the results of a numerical simulation, to efficiently perform the inverse analysis. Even though performing MC on a metamodel of the process is fast, the uncertainty propagation step can be performed more efficiently using an analytical approach [11]. Therefore, in this work analytical propagation of uncertainty is combined with metamodeling to perform the inverse analysis more efficiently.
To the best of the authors' knowledge, this kind of problem has not been addressed in the literature and more specifically in the forming-related literature. In this article, a meta-modelbased approach is implemented and the goal is to specify the output tolerance and tailor the noise scatter. This article is structured as follows: In Tailoring Input Noise Scatter, the robust optimization approach and tailoring noise scatter is introduced and formulated. In Application in Metal Forming Process, the numerical simulation of a lab-scale B-pillar is introduced and the input and output of the process for both the forward and the inverse method are defined. The results are presented and discussed in Section Results and Discussion. Figure 1a shows the concept of robust optimization schematically. In this figure, x represents a design parameter, z represents a noise parameter that has a known probability distribution, and r is the response. The goal is to find an optimum design, x opt in which the variation of response around the target mean is minimized.

Forward problem
For this purpose, the relation between r and (x, z) is expressed using a metamodel. To obtain the metamodel of the process, a design of experiments (DOE) on the (x, z) domain is generated and r is evaluated at each DOE point. The evaluation is usually regarded as a result of a blackbox function evaluation. For optimization of metal forming processes, a finite element (FE) model can be used as a black-box function.
A metamodel is then fitted to the responses obtained from this black-box function evaluation [14]. In the robust optimization method the search for an optimum design, x opt , is performed based on propagation of the given noise through the metamodel. An objective function is then defined in terms of statistical parameters of response to evaluate the robust optimum.
In a general case, the vector of design variables, x, and the vector of noise variables, z, are used to evaluate the response of the process, r(x, z). The objective function can be defined in many ways such as the variation of response (σ r (x)) the deviation of mean of the response (μ r (x)) from a target value (C r ) or a combination of these two [5]. The latter is used in this article and the objective function aims to minimize the variation of the response while setting the mean value of the response on target value: where w is a weighting factor to adjust the optimization objective between mean on target and output variation. μ r c Fig. 1 (a) Finding robust optimum based on uncertainty propagation from a given input noise (b) Finding acceptable noise from specified output tolerance is the mean of the response of constraint and σ r c is its standard deviation. The choice of n for constraints is related to the desired reliability level. It is generally assumed that the noise input follows a known normal distribution that has a fixed or nominal mean and standard deviation [15]. Therefore, any change in the input mean and standard deviation leads to a different solution for the robust optimization. For this reason, Eq. 1 is rewritten by: Eq. 2 reflects the fact that the objective function is a function of the input mean and standard deviation while in robust optimization it is evaluated using nominal values.

Inverse problem
The concept of tailoring input scatter is illustrated in Fig. 1b. If the output variation at the robust optimum obtained using robust optimization does not meet the specified tolerance, the input noise must be adjusted. In this approach, one works back from the specified response variation to the acceptable input noise. Considering the non-linear behavior of the process, the presence of multiple noise parameters and the interactions among them, an efficient approach must be developed. Ideally, the input noise must be zero to achieve zero variation in the response. However, a tighter control of noise parameter incurs additional costs. Therefore, to minimize costs, thevariation of the input noise parameters can be accepted up to a certain level where the response meets the specified tolerance. The inverse approach for tailoring the noise scatter is formulated as: whereσ z are defined as σ z /σ nom z . The objective is to find a set of (σ z , x) that maximizes the weighted sum of the normalized noise variation. Besides, the constraints are assessing the feasibility of the obtained solution. A weighted sum formulation is used as it is often employed to provide feasible solutions by varying the weighting factors [16]. The objective function in Eq. 3 handles the following requirements simultaneously: wider specifications for noise parameters reduce the costs of the process. Therefore, maximizing the sum of noise variations is required to obtain a cheaper process. Various noise parameters can have different dimensions, and therefore, each noise parameter is normalized using its nominal value. The weighting factors are associated with the relative cost of control for each noise parameter. It means that if a noise parameter is more costly to control than other parameters, it gains a bigger weight; thus the maximization problem concentrates more on maximizing the variation of that noise parameter. Based on this explanation, the weights are non-dimensional coefficients that are obtained using the costs of controlling each noise parameter. As an example, assume that z 1 and z 2 are two noise parameters. If reducing σ z 1 and σ z 2 by 50% cost 50A C and 100A C, respectively, the weight factors are ω 1 ∝ 50A C and ω 2 ∝ 100A C. When they are normalized such that ω i = 1, ω 1 = 1/3 and ω 2 = 2/3 are obtained.
It must be noted that the objective function in Eq. 3 can be formulated using other expressions, and the choice of objective function is not limited to the above-mentioned objective function. In Eq. 3, the constraints handle the feasibility of satisfying the specified tolerance. It can be seen that the objective of the robust optimization in Eq. 2 becomes the constraint of the inverse problem in Eq. 3. In addition, it is assumed that the noise mean is fixed on the nominal value and the focus is on the effect of the variation of noise parameters.
Compared to Eq. 1 in which the noise mean and standard deviation are fixed and the search for robust optimum is performed only on x, the proposed method is computationally more expensive. Even though a metamodel of the forming processes is employed to perform the inverse analysis, the computational effort is significant. In this article, instead of conventional techniques such as Monte Carlo (MC) and its variations, an already developed analytical method for robust optimization [11] is employed to calculate the uncertainty and optimum in each step of inverse analysis. It has been confirmed that using the analytical method during robust optimization procedure improves the speed and accuracy of the calculations compared to existing approaches.
Using the results of univariate integrals to evaluate multivariate integrals and obtain the output moments [17] the mean and standard deviation of the output are calculated using: In these equations, a 0 is a constant, a i and a j are the weights obtained by fitting a metamodel, n x is the number of design parameters, n z is the number of noise parameters, N is the number of DOE points, and b ip (x p ) and b jp (x p ) are the basis functions of the metamodel. These basis functions are defined by: where x p is the evaluation point, x ip are the DOE points and the values for ϑ p are obtained from fitting the Kriging metamodel. In this approach, the data is assumed to be a realization of a Gaussian random field and therefore spatially correlated. The correlation between the poitns depend on their distance from eachother. A normal probability distribution is defined by: where μ q and σ q are the mean and standard deviation of the input noise data. Then C1 iq and C2 ij q are obtained for a Kriging metamodel and a normal probability distribution using:

Finite element simulations
In car body engineering, the pillars between windows and/or doors are designated alphabetically from the front: A-, B-C-pillar. The B-pillar is between front and rear doors and plays an important role in side impact crash safety. It needs to have a reinforcement that is strong and deep, putting high demands on formability. In the lower part of the Bpillar the forming problems are largest so this geometry can be considered a critical test case. Hence, this (T-shaped) geometry was selected for the optimization problem. FE simulations have been validated in a previous study by Abspoel et al. [18] and the FE results have shown a good agreement with the experiments.
In the lab-scale B-pillar forming process, a metal sheet is stretched and bent over the die using a punch. The metal sheet made of dual-phase (DP 800) steel is placed on the die. The blank holder moves towards the metal sheet to clamp it in position. Then the punch stroke shapes the sheet. The FE model is built in AutoForm R7 software and its simulation results are used to make a metamodel of the process. The configuration of the FE model is shown in Fig. 2. Elasto-Plastic Shell elements using 11 integration points in the thickness direction (EPS-11) are used for the simulations. The design and noise parameters of this process are specified in Table 1. The design parameters which can be adjusted during forming are blank-holder force (F bh ) and blank corner radius (R). According to this table these The design parameters are shown in Fig. 2 in which F bh is the force applied uniformly on the metal sheet.
The material behavior was obtained using a tensile test combined with a stacked compression test. For the yield locus, Vegter 2017 was used as described in [19]. Uniaxial stress factors and Lankford coefficients are needed to generate an accurate yield locus using this method. Those parameters can be determined only form tensile tests, and therefore Vegter 2017 is suitable for variation studies. Swift model is used to describe the hardening of the material by: where σ is the stress, ε is strain, and K, n, and ε 0 are the parameters that are usually obtained from uniaxial tensile tests. In this work, ε 0 = 0.002 and n = 0.169. In addition, the Coulomb friction model is used in the simulations. Three noise parameters are identified in this process which are the strength coefficient in the Swift hardening law (K), the friction coefficient (m) in Coulomb friection model, and the thickness of the sheet (t). It is assumed that the noise parameters have a normal probability distribution and that there is no correlation between those parameters. To use Equations 4 and 5, it is assumed that the noise variables are statistically independent. In the current approach, K depends on the processing history and the elemental composition of the material, t is usually controlled in the last rolling step during sheet production stage, and m originates from surface texture and lubrication condition during the process. The deformation in the last rolling step is so small that the influence of thickness reduction on hardening behavior can be ignored. Therefore, it is a rational assumption that the three input parameters can vary independently of each other. After releasing the force and removing the part from the die, a representative final shape of the part is shown in Fig. 3. The main response is the final angle of the profile (θ ), and a constraint is defined through the percentage thinning (δ) in the corner region. Preliminary simulations confirmed that the order of variation due to numerical errors in FE simulations is lower than the order of variation caused by variation of noise variables. This step is recommended by Wiebenga et al. [20] to ensure that the results obtained from robust optimization are valid.

Sampling and constructing a metamodel
The Latin hypercube sampling (LHS) technique is used to create a DOE consisting of 90 points in the five-dimensional parameter domain. The responses, θ and δ, are calculated using the FE simulations. Figure 4 shows the final shapes of the profiles and the variation of θ , on DOE points obtained using FE simulations. The Kriging method is used for metamodeling as it is commonly used in forming processes to model the non-linear relationship between input and output [21]. It is used during both robust optimization and tailoring the material scatter. The analytical method [11] is used to calculate the mean and standard deviation of the responses during analysis.

The forward problem
Based on Eq. 1, the forward problem is defined as: Throughout this chapter, the minimum objective function of forward problem is referred to as the process quality indicator. A forward analysis is performed using nominal set of noise parameters to find the nominal process quality indicator which is denoted by C nom .

The inverse problem
In the inverse method, the required process quality indicator is chosen and the goal is to find the maximum variation in the input noise that satisfies the requirements of the output. The inverse approach is expressed by: Using this formulation, different requirements on response can be set using parameter C tol , and a single objective function is defined using weighting factors, ω i , for each noise parameter. The presence of the upper bounds for σ K , σ t , and σ m in Eq. 12 are due to the extension of the DOE on which the metamodel is based. When one of those becomes active, it means that the DOE must be extended in that direction. This is discussed in Section Intact Output Tolerance, Wider Specifications.

Results and Discussion
The result of robust optimization using nominal noise variations is obtained based on the analytical approach presented in Tailoring Input Noise Scatter and the results are presented in Table 2. At the robust optimum design, the  variation of θ is shown using a normal distribution in Fig. 5. This figure represents the minimum scatter at the optimum design setting, in which the weighted sum of 'mean on target' and 'response variation' is minimum. The process quality indicator is obtained at that specific design setting and is presented in Table 2. For comparative purposes, the results of the MC sampling method are also shown in Fig. 5 in a bar plot. MC result shows that the main response closely follows a normal distribution and the mean and standard deviation of the output obtained via analytical calculations can be reliably used to accurately predict the robust optimum. A slight deviation from normal distribution is observed in this figure, however, it is negligible. For a large deviation from a normal distribution, one can employ higher order moments of the output during robust optimization to evaluate the robust optimum more accurately [22]. For tailoring the scatter two scenarios are defined. First, the process quality indicator is set below the nominal process quality indicator, C tol < C nom , and the inverse problem is solved to find the acceptable variation of noise variables which will be less than nominal values. In the second scenario, the process quality indicator remains intact, C tol = C nom , and the inverse problem is solved to relax the variation of noise variables. These scenarios are addressed in the following sections.

Tighter Output Tolerance, Tighter Specifications
As a starting point, one can set the required process quality indicator on C tol < C nom and solve the inverse problem. The results of an inverse analysis for various values of C tol are shown in Fig. 6. In all cases, the weighting factors are equal to [1/3, 1/3, 1/3], meaning all noise factors are considered equally expensive to control. This figure shows how reducing the value of C tol , decreases the acceptable variation of each parameter. For the majority of the range of C tol , the required standard deviation of the friction coefficient remains at its specified upper limit. This means, that a relaxation of this upper limit would give a cheaper optimum. In addition, the acceptable variation of each parameter does not proportionally decrease by reducing C tol due to the non-linear response to the variation of the parameters and different sensitivities. Figure 7 shows the acceptable variation of noise parameters for C tol = 0.1C nom and various weighting factors. This figure shows that if a parameter is more costly to control than the other ones, the accepted variation for that parameter increases, however, to achieve the specified tolerance another parameter should be controlled tighter. Table 3 shows the optimum settings in these cases. This table shows that, depending on the settings in the inverse approaches, different process design parameters are obtained than when the forward approach is used. In some cases the bound constraint on R is active, and therefore the solution of the optimization problem is limited by this bound. Nevertheless, the bound constraint on F bh is never active. Outside the bounds for design parameters, the metamodel is not accurate. In these cases the metamodel can be extended in that specific direction. In the next section, it is shown how to extend the metamodel based on a solution obtained from the inverse problem.
Pareto optimal solution sets can be obtained for multiobjective optimization problems to show the trade-off between the objectives [23]. In this work only one objective function is used comprising the weighted sum of the variation of various noise parameters. Therefore, a variety of optimal solutions can be found by modifying weighting factors to meet the constraints. Fig. 8, shows the optimum acceptable variation of noise parameters when C tol = 0.5C nom . In this case, it is not required to controlσ m and it varies within the nominal range. Therefore, the Pareto front is plotted only in 2D. Despite the weighting factors being selected at euqal intervals, the points on the Pareto front are  [24]. Since the shape of the Pareto front is not known in advance, it is not feasible to determine the values of ω that lead to a uniform distribution of the points on the Pareto front.
As mentioned earlier, when C tol is reduced to 0.1C nom , all three noise parameters must be controlled tighter. For that reason, a 3D Pareto front is plotted and shown in Fig. 9. This surface represents the fit over 55 solution sets and demonstrates the trade-off among the variation of noise parameters.

Intact Output Tolerance, Wider Specifications
In the previous section, C tol was set lower than the minimum objective function in robust optimum, C nom , meaning that a tighter tolerance of the output is required. In this section, Table 3 The acceptable variation of noise and the optimum design for various weighting factors  9 Fitted Pareto front for C tol = 0.1C nom by setting C tol = C nom the aim is to achieve a cheaper process specification that results in the same process quality indicator. This usually means that, σ exceeds the bounds in Eq. 12. Theoretically, the noise variation can be increased outside of those given bounds since a material or process with wider specifications is cheaper. However, these bounds are present to avoid using the metamodel of the metal forming process beyond those bounds. As there are no DOE points outside the ±3σ nom bounds, the metamodel of the process is less accurate in those regions. Therefore, those regions are omitted to avoid inaccurate predictions occurring by extrapolating the results. Nevertheless, the possible solution is to add new DOE points to improve the prediction of the metamodel in those regions.
First of all, it is required to realize which noise bounds will be violated in the case of C tol = C nom . For this purpose, the initial metamodel, despite being not accurate outside the bounds, is used to estimate the required extension of noise bounds. An inverse analysis using equal weights for all noise parameters considering Eq. 12 is performed without including the bounds for σ . The results of the preliminary analysis are shown in the first row of Table 4. As expressed in this table,σ K andσ t are within the initial noise variation, butσ m exceeds the initial limit. Since there are no DOE points from 3σ m to 3.8 × 3σ m and from −3.8 × 3σ m to −3σ m , the results of the robust optimization and inverse analysis is poor. As observed in other studies [25], a more complete sampling, around the periphery and within the study region is required to make an accurate metamodel. Therefore, new DOE points are added to the regions without sampling points.
The strategy of adding new DOE points is explained using Fig. 10. An initial 2D DOE, in which z 1 and z 2 represent two noise variables, is shown schematically in Fig. 10a. Assuming that an inverse approach using the initial metamodel finds a solution in which σ z 2 exceeds σ nom z 2 (σ z2 > 1), the DOE can then be extended using infill points in that direction. A larger variation of one noise will lead to a smaller variation of another noise. This has been shown in Fig. 10b using the dashed lines. Therefore, the new infill points are confined in the range of ±3σ z 1 in z 1 direction. The new infill points are shown by circles in Fig. 10b.
According to the first row of Table 4, the DOE should be extended in m variable direction. The initial DOE is enriched by adding six infill points. The LHS method is Table 4 The acceptable noise variation and the optimum design outside the initial noise limits using improved metamodel  Fig. 10 The strategy of adding infill points used to generate 5000 sets of the new infill points in the expanded region, from which a set having the maximum sum of the distances from initial DOE points is selected. The FE simulation is then performed on those newly added points and the procedure of inverse analysis is repeated using the metamodel built on the updated DOE. Table 4 shows the results of using the updated metamodel in the inverse approach. In this table, adding only six infill points improves the results of the inverse analysis significantly. In addition, the infill points depend on weighting factors used in the analysis. Therefore, a change in the weighting factors might lead to the change of the direction of recommended expansion of the metamodel.

Outlook
In this work, a general method for tailoring of scatter is presented and implemented. Although specific assumptions have been made for a case of metal forming process, this method is flexible and it can be modified to be applied in other engineering problems.
First, the objective function is chosen to maximize the weighted sum of normalized sigma. In this work, the objective function is normalized using the nominal variation of noise. Using this method, the percentage acceptable variation is obtained for each parameter during tailoring material scatter. However, the normalization can be performed using various methods [26]. Nevertheless, normalization is an important step to make it feasible to set the weights such that they are significant relative to each other and relative to the objectives [27]. It is also conceivable that the user defined input of the cost parameters is in a readable form (e.g. A C per mm thickness control) and the normalization is done inside the software during the analysis.
Second, in Eq. 3, the material scatter is tailored considering a fixed mean value for all noise parameters.
When the mean value for some applications needs to be adjusted, a new objective function must be defined to take into account the costs of the deviation of the mean value of each noise parameter from the nominal mean. Therefore, the objective function must be able to maximize not only the noise scatter, but also the deviation of noise mean from the nominal mean value.
Third, the first inequality constraint can be considered as an equality constraint meaning that it is always active within the computational accuracy. This is due to the fact that inactivity of this constraint means that there is still room for increasing the variation of the noise parameters. It is not expected that a material with tighter tolerance could be cheaper to produce than a material with larger variations.
Last but not least, the objective function in the inverse method is a linear combination of noise variations to be maximized. Determining a set of weights for the relative importance of objectives is only locally valid [27] and they cannot reflect the relative importance of objectives in the whole domain. For instance, reducing the scatter of a noise parameter by 50% can be relatively cheap, whereas, it might be very costly to reduce it by 25% due to nonlinear nature of controlling a noise parameter. Therefore, the weights and their applicable domain must be selected carefully or a nonlinear weighting factor must be selected.

Conclusions
An inverse robust optimization approach is presented and implemented to tailor the variation of material and process noise parameters based on the specified product tolerance. This method has been successfully applied to the forming process of a lab-scale B-pillar. The results show how the acceptable material scatter varies by setting a tighter product tolerance or setting various weighting factors. This approach is flexible in terms of adjusting weighting factors which are related to the costs of controlling each noise parameter. Not all noise parameters need to be controlled up to a certain reduction of the product tolerance. However, further reduction of the product tolerance, requires more control on all noise parameters. For these cases, the Pareto fronts are presented to show the influence of employing various weighting factors and demonstrate the trade-off between controlling the parameters to minimize costs. Additionally, the developed method could equally be used to relax material specifications when a process meets the initial tolerance by a large margin. This, however, requires a re-construction of the meta-model on a wider (noise) base to ensure that it is sufficiently accurate.