Using Bayesian optimization for warpage compensation in injection molding

In injection molding, shrinkage and warpage lead to a deformation of the produced parts with respect to the cavity shape. One method to mitigate this effect is to adapt the cavity shape to the expected deformation. This deformation can be determined using appropriate simulation models, which then also serve as a basis for determining the optimal cavity shape. Shape optimization usually requires a sequence of forward simulations, which can be computationally expensive. To reduce this computational cost, we use Bayesian optimization which uses Gaussian process regression as a reduced order model. Additionally, Gaussian process regression has the benefit that it allows to account for uncertainty in the model parameters and thus provides a means to investigate their influence on the optimization result. We present a Gaussian process regression trained with samples from a finite‐element solid‐body model. It predicts the deformation of the product after solidification and, together with Bayesian optimization, allows for efficient cavity optimization.


| INTRODUCTION
Injection molding is a widely-used manufacturing process for plastics products.In order to ensure a high part quality, a number of causes for inaccuracies need to be considered.One main cause for reduced shape accuracy is warpage; a result of both uneven shrinkage and an inhomogeneous temperature distribution within the part during and after solidification [1].
Warpage is influenced by a variety of factors and numerous studies have investigated methods to reduce warpage by adjusting process parameters or part design [2].In that work, the influence of process parameters such as melt temperature, injection time, injection This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.pressure, holding pressure, and cooling time is explained and different optimization strategies are discussed.The optimization techniques include genetic algorithms, response surface methodology, Bayesian optimization, and particle swarm optimization.In terms of component design, the wall thickness of the part can be optimized to reduce warpage [3].Also, the design of the cooling system can reduce the warpage [4].
However, these methods cannot completely eliminate warpage.A potential approach to further minimize warpage is to compensate for it by adjusting the cavity shape.In this way, the part is produced with a deformation that then compensates for the warpage that occurs.
Finding the optimal cavity shape that compensates for warpage is difficult.The simplest method is based on experiments and trial runs.The cavity is iteratively modified based on measured data from parts produced in trial runs.For this method, it is critical to use enough measurement points with sufficient accuracy.It is also possible to compute the corrected cavity shape based on measurement data of the entire surface of the part [5].
The industry also uses simulation-based compensation methods.With commercially available simulation software, the warpage can be predicted.Integrated optimizers and sensitivity analysis for process parameters and simple geometric features are also available [6,7].
An automated method to compute the optimal cavity shape is still a work in progress.Several methods have been proposed in the literature.
A method called 3D volume shrinkage method has been used for shrinkage compensation by scaling the part in each direction of the coordinate axis, resulting in a reduction of the warpage [8].
Another method is to use an inverse model to compute the optimal cavity geometry [9,10].A forward simulation is run once to obtain the temperature distribution in the part and then an inverse model is employed to identify a modified part that after shrinkage and warpage corresponds to the desired product.This geometry is then used as the new cavity shape.The advantage of this method is that only one run of the inverse model is required.The disadvantage is that the temperature distribution is computed only once and is not adjusted for the compensated geometry.
In another iterative method, the objective function is the difference to the desired geometry and it is computed locally [11].With this local objective function, the geometry can be adjusted locally to reduce each objective function with each iteration.This approach requires only a few iterations.Another advantage of this method is that any forward simulation model can be used.
In this paper, we propose a shape-optimization-based approach to finding the cavity shape that minimizes warpage.In general, shape-optimization approaches require four ingredients: (1) a shape parameterization, (2) an easy-to-evaluate forward simulation, (3) an objective function, and (4) an optimization algorithm.In terms of shape parameterization, a low-dimensional representation has the advantage that the computational effort is more manageable and therefore the set of suitable optimization algorithms is larger.Furthermore, such parameterization may facilitate the compatibility with computer-aided design systems.There are several ways to achieve such a parametrization [12].Here, we use the so-called Free-form deformation method.Free-form deformation not only allows for local adjustments of the geometry, but it also has the advantage that it is easy to set up also for complex geometries as it is based on the embedding of point clouds in a box spline, which we emphasize is not boundary-conforming to the represented geometry [13].
The geometry can then be modified by displacing the control points of the box spline.The choice of spline type for the box spline can significantly affect the flexibility and computational efficiency of the optimization process.Common spline types include B-splines, nonuniform rational B-splines, and Bézier splines [14,15].In our study, we chose B-splines because they offer a suitable balance between flexibility and computational cost, providing a smooth and continuous representation of the cavity geometry while maintaining computational efficiency [14].
Since shape-optimization problems require many queries of the forward simulation, as a second ingredient we require an efficient forward simulation [16].In this work, we will use Bayesian optimization, a surrogate model-based optimization algorithm.Making use of surrogate or reduced-order models can serve as an important means to lower computational complexity [17].Surrogate models and reduced order models are especially useful in optimization problems where the objective function is expensive to evaluate [16].This is the case, for example, when the objective function is computed by running a full-order finite element simulation model.
In this paper, we generate a surrogate model that represents the functional relationship between the optimization parameters i. e., the free-form deformation control point positions and the objective function.
Such surrogate models can be obtained in many ways, including regression techniques or machine learning [16][17][18][19].In contrast to reduced-order models, surrogate models are purely data-based.This means that they are non-intrusive and furthermore can be trained in a hybrid way using both simulation and experimental data [16].In this work, we obtain the surrogate model using Gaussian process regression.The reason for this choice is that in addition to the surrogate model itself, Gaussian process regression also provides an uncertainty bound for its prediction of the objective function [20].Knowing how accurate the Gaussian process regression model is for certain optimization parameters can be used during training to calculate the optimal next training point.This is often referred to as Bayesian optimization.To find the optimal next training point, an acquisition function is defined that balances exploitation and exploration, i. e., the weighting of low function value versus high uncertainty.Thus, promising points with low function value and high uncertainty are chosen, as here the probability of finding a better point is highest [21,22].Note that in literature, Bayesian optimization has been applied to finding optimal process parameters for injection molding [23][24][25].The difference to this work is that there the process parameters such as the cooling time or the holding pressure were optimized, whereas in our work, the cavity shape is optimized to compensate for warpage.
In classic Bayesian optimization, the bounds for the optimization variables must be set first.If no bounds are set, the learning function would propose training points that are too far apart in the search space, which can be a drawback if the region where the global optimum lies is unknown.Several methods have been proposed in the literature to deal with this problem.One method is to constrain the uncertainty during optimization with an upper bound [26].Another method is to use initial bounds and expand them when the bound is reached, called adaptive expansion [27].
The paper is organized as follows.In Section 2, we describe our proposed method.This is followed by a numerical example in Section 3 and the discussion and outlook in Section 4.

| A BAYESIAN OPTIMIZATION FRAMEWORK FOR IDENTIFYING THE OPTIMAL CAVITY SHAPE
In this section, we detail the methods used in our study, including the spline-based free-form deformation and the simulation model, followed by the objective function and the implementation of Bayesian optimization to ultimately achieve the goal of warpage minimization by finding the optimal cavity geometry.
In this paper, surrogate-model-based optimization is used to solve the geometry compensation problem.Other components involved are the objective function, which is computed as the mean distance to the desired geometry.The adaptation of the cavity geometry during the optimization iterations is achieved via a spline-based free-form deformation.

| Spline-based free-form deformation
To adjust the cavity geometry in a flexible and efficient way, we use a spline-based free-form deformation approach [13,15,28].In particular, our deformation spline is a B-spline.The deformation spline is usually a boxshaped spline, i. e., either a rectangle in 2D or a cuboid in 3D.In many cases, the control points are homogeneously distributed, but this is not a strict requirement.The finite element mesh representing the initial geometry, i. e., the cavity geometry in the warpage simulation is then embedded in the box spline.In simple terms, this means that each finite element node is attached to a specific point on the spline.As the spline is deformed during the optimization iterations, also the finite element mesh deforms accordingly and the cavity shape is modified.In practice, this means that for each finite element node one identifies a local coordinate on the spline.When the spline is re-evaluated after its deformation, this automatically identifies the new positions of the finite element nodes.This means that no remeshing of the domain is necessary.By manipulating the positions of these control points, the cavity geometry can be deformed in a smooth and continuous manner.
It is to be expected that our approach requires more iterations than the local objective function method, since the optimization problem is of higher dimensions.In particular, the dimension is computed as the number of spline points multiplied by their respective number of degrees of freedom.Depending on the number of spline points used, the dimension may increase rapidly.However, given the high both global and local shape flexibility offered by free-form deformation along with the guarantee of smooth shapes, we consider this additional computational effort justified.

| Simulation model
The injection molding process can be subdivided into several phases, namely, filling, packing, solidification, and demolding.Our warpage model sets in at the process stage when the component is fully solidified, but not yet cooled down to room temperature [29].This is the main phase where warpage occurs as the component is no longer constricted by the mold.As initial data, the model requires the temperature distribution within the component after the solidification.The inhomogeneous temperature distribution is the driver of warpage and can be determined using a separate simulation that is not part of this paper.
As constitutive equation, we use a steady, linear thermoelasticity model [30].The domain where we solve our model is the volume of an injection molding part, denoted by Ω.As discretization, we use a tetrahedral mesh with 29048 elements and 8725 nodes.To prevent translational or rotational movement, in total six degrees of freedom near the gate are restricted as the boundary conditions.The stress tensor is computed as: where C is the stiffness tensor and ɛ the strain tensor.These are defined as follows with the identity matrix I and the fourth order identity tensor J: The stiffness tensor requires the Lamé parameters λ and μ, which can be defined in terms of Young's modulus E and Poisson's ratioν.
Additionally, u is the deformation, α is the coefficient of thermal expansion and T is the temperature for the initial and final state respectively.
The model then calculates the steady-state solution of the deformation when the entire part has reached ambient temperature.For this, the weak form with test functions ϕ k of the momentum balance is solved over the domain Ω.
The simulation model thus provides a prediction of the distortion that will occur after the cooling phase.The deformed geometry is then used in order to compute the objective function for the optimization process.It is evident that this model makes a number of assumptions, such as linearity, unrestricted cooling, etc.However, the correctness and accuracy of this simulation model is not the focus of this study.In general, any warpage prediction model can be used.

| Objective function
For optimization, the calculated warpage must be quantified by an objective function.The objective function measures the deviation between the warped geometry and the ideal geometry of the molded part.In our study, we define the objective function as the average distance between corresponding points on the warped and ideal geometries.
We use this point-based distance measurement, as this is a very simple method which is easy to implement and sufficiently accurate for our use case.The points are sampled uniformly over the surface of the undeformed part, based on the mesh elements.As the same mesh connectivity is used for the deformed and undeformed mesh, each point on the deformed domain has a corresponding point on the undeformed domain.We calculate the Euclidean distance between each pair of corresponding points and sum these distances divided by the total number of points: with x i 2 G � @W and xi ðsÞ 2 GðsÞ � @ WðsÞ Here, G or G is the set containing the sample points and @W or @ W denotes the surface of the undeformed or deformed part, respectively.Note that the shape of the deformed surface depends on the coordinates of the spline control points s.The coordinates of the part are x for the ideal geometry and x for the deformed geometry.This distance-based metric provides a comprehensive measure of the overall distortion in the part and allows direct comparison of different cavity geometries during the optimization process.With this objective function, the optimization problem can be formulated as follows: where s is a flattened vector representing the coordinates of the control points, combining the values of the coordinates from each direction.The lower and upper bounds of the spline control point coordinates are given by the vectors B lower and B upper .The condition between the vectors s, B lower and B upper has to be read componentwise.By minimizing this objective function, we aim to reduce the overall distortion in the molded part, thereby improving its quality and adherence to the desired specifications.

| Bayesian optimization implementation
Using Bayesian optimization, we identify the optimal cavity geometry in the sense of Equation ( 6).The Bayesian optimization is based on a Gaussian process regression surrogate model determined using the opensource Python library psimpy.For the Gaussian process regression, we use a Matérn kernel with the parameter ν Matérn = 2.5.The expected improvement function EI is used as the acquisition function for the Bayesian optimization to determine the most promising next candidate geometry in the search space [31].The expected improvement function is defined as: with Here, m is the mean and σ is the variance of the Gaussian process regression prediction.F(s + ) is the best observed objective function value, Φ(Z) is the cumulative distribution function of the standard normal distribution, and ϕ(Z) is the corresponding probability density function.The parameter ξ balances exploration and exploitation.
The best candidate geometry lies at the minimum of this exploration function.Thus, in each iteration step, the (ideally global) optimum of this function is found.This optimization problem is: For an efficient Bayesian optimization, especially in high dimensions, it is important not to have a too large search space for the global optimum.Otherwise, too many iterations will be needed for the Bayesian optimization to converge.Therefore, we start with small bounds for the optimization variables, which are expanded adaptively.The criterion for expansion is that the current optimum lies at the current boundary.Then the boundary is expanded in that direction.With this method alone, the algorithm tends to over-explore the domain, leading to unnecessary simulations in regions where the objective function is far from the optimum.
Therefore, in addition to the adaptive expansion, the search space is limited by imposing an upper bound on the optimization of the acquisition function.This results in the adaptive equation of the acquisition function optimization: Here, τ is an additional parameter that reduces the exploration of the Bayesian optimization algorithm.To find the global optimum of this constrained problem, the Python library scipy.optimize is used.The Basin-hopping algorithm is used, which starts several local optimizers in the search space.The local optimizer uses the sequential least squares programming algorithm.For both optimizers and the Bayesian optimization algorithm, multiple hyperparameters, as mentioned in this section, are selected, Table 1.Here, η and the bounded expansion rate have the most influence on the convergence of the algorithm.The result of this optimization is a set of spline control point coordinates.
The Bayesian optimization algorithm starts with the initial cavity shape and the warping model is used to compute the objective function, Figure 1.The objective function provides the data to train the Gaussian process regression.The Gaussian process regression is used to find the minimum of the acquisition function.The location of the minimum of the acquisition function selects the coordinates of the spline points.Free-form deformation deforms the mesh of the next cavity shape, and the loop starts again.The algorithm then runs for a fixed number of iterations.

| NUMERICAL EXAMPLE
In this section, we present a numerical example for our optimization approach, highlighting the effectiveness in minimizing warpage in injection molded parts.As stated in the Section 2, the warpage model requires an initial temperature distribution.This distribution was determined using the commercial software Moldflow.For this simulation multiple process parameters need to be selected, Table 2.In the subsequent warpage simulation according to Section 2.2, the necessary material parameters were chosen, Table 3.
The geometry is a slice of a part with multiple fins and the initial temperature field is indicated by the color field, Figure 2. As such, it can be reduced to a 2D-geometry.For symmetry reasons, only half of the part is considered.The temperature just after solidification is unevenly distributed.The highest temperature is at the injection point at about 520 K and drops rapidly for the rest of the part.The tips of the fins feature the lowest temperature of 360 K.This image also shows the spline control points used for the free-form deformation.To keep the computational cost low, we use as few control points as possible.The underlying B-spline is comprised of nine control points, whose positions are all subject to optimization.This results in a total of 18 degrees of freedom, namely the coordinates of the spline control points in x-and y-direction.These spline control points are then moved with a sample distribution to deform the geometry, Figure 3.
The optimization algorithm was run for 300 iterations using the desired geometry as the initial cavity shape.The objective function value indicating the quality of the cavity shape decreases with increasing number of iterations, Figure 4.The initial value of the objective function is about 2.0 mm.In the first 50 iterations, the objective function value decreases rapidly to about 0.5 mm.In the next 250 iterations, the objective function value decreases much slower, but distinct improvements can still be observed.For comparison, also the SciPy implementation of the sequential least squares programming (SLSQP) algorithm was run for the same problem.The decrease of the objective function value is slightly worse compared to the Bayesian optimization algorithm.
In the case of our simple numerical example, the sequential least squares programming (SLSQP) algorithm did not get stuck in a local minimum.However, in a more complex case, local minima may be a problem, in which case Bayesian optimization is more likely to find the global optimum.The optimization algorithm results in an improved cavity shape, Figure 5.By choice, the initial cavity shape corresponds exactly to the desired product geometry, Figure 5a.As such, this geometry has a dual-purpose.When the initial cavity is used, the right corner of the part is bent up and the two largest ribs are bent slightly inward, Figure 5b.For the optimal shape determined via Bayesian optimization, the right corner of the cavity shape is moved up and the ribs are bent inward, Figure 5c.Notice, that roughly the opposite deformation of the deformed part is applied.The final part resulting from this cavity shape shows that the deformation has been significantly reduced, Figure 5d.Some minor deviations compared to the desired product geometry remain, for example the length of the third fin is too short.

| DISCUSSION AND OUTLOOK
As demonstrated by means of the numerical example, the introduced algorithm is able to improve the cavity geometry in such a way that the warpage is significantly reduced.The key advantage of our method is that it is non-intrusive, meaning that any warpage model or experimental data can be used to generate the surrogate model and optimize the cavity shape.In addition, this method is easily adaptable to include additional optimization parameters for warpage minimization, such as process parameters.
However, in comparison to the local objective function method it should be noted that our method requires more than an order of magnitude more iterations.This is due to the fact that our method uses a global objective function that averages the distance to the desired geometry over the entire part.With this method, the information about where the deformations occur on the part is lost and must be "rediscovered" by acquiring enough training data from the simulation model.It is necessary to find out which spline point corresponds to an area with larger deformations.Depending on the number of spline points, this connection is difficult to find.This will be addressed in future work.
Furthermore, we currently do not restrict the possible range of cavity shapes except for the imposed bound constraints on the displacement of the spline control points.As a result, the adapted cavity shape may result in a produced part that is difficult to demold, i. e., difficult to extract from the cavity.While curved parts are not impossible to demold, it is more expensive.As a remedy, one could implement an appropriate constraint on the cavity shape.It is also important to note that no process parameter optimization was performed in the example shown.With optimized process parameters, the total warpage will be reduced and thus a demolding criterion may no longer be necessary.

1 F I G U R E 1
Overview of the Bayesian optimization algorithm.

T A B L E 2 À 1 F I G U R E 2 F I G U R E 3
Injection molding parameters for Pocan B 1305.Geometry of the part, including the temperature field and the spline control points.Deformed geometry (gray) with the corresponding spline surface and control points (black).

F I G U R E 4
Convergence of the Bayesian optimization algorithm and the sequential least squares programming (SLSQP) algorithm.On the y-axis, the objective function value is shown, which measures the average distance of the deformed part to the desired part.On the x-axis the number of iterations is shown.F I G U R E 5 Comparison of the initial cavity shape (a) and the optimized cavity shape (c) with their corresponding deformed parts (b) and (d).