Multi-Objective Optimization of Plastics Thermoforming

: The practical application of a multi-objective optimization strategy based on evolutionary algorithms was proposed to optimize the plastics thermoforming process. For that purpose, in this work, differently from the other works proposed in the literature, the shaping step was considered individually with the aim of optimizing the thickness distribution of the ﬁnal part originated from sheets characterized by different thickness proﬁles, such as constant thickness, spline thickness variation in one direction and concentric thickness variation in two directions, while maintaining the temperature constant. As far we know, this is the ﬁrst work where such a type of approach is proposed. A multi-objective optimization strategy based on Evolutionary Algorithms was applied to the determination of the ﬁnal part thickness distribution with the aim of demonstrating the validity of the methodology proposed. The results obtained considering three different theoretical initial sheet shapes indicate clearly that the methodology proposed is valid, as it provides solutions with physical meaning and with great potential to be applied in real practice. The different thickness proﬁles obtained for the optimal Pareto solutions show, in all cases, that that the different proﬁles along the front are related to the objectives considered. Also, there is a clear improvement in the successive generations of the evolutionary algorithm.


Introduction
Thermoforming is a thermoplastic processing technique commonly used in the rigid packaging industry. A variety of thermoplastic materials can be used in this process, including semi-crystalline polymers, such as High Density Polyethylene (HDPE) and Polypropylene (PP) and amorphous polymers, such as Acrylonitrile Butadiene Styrene (ABS), Polystyrene (PS) and High Impact Polystyrene (HIPS). HIPS is a lightweight and inexpensive thermoplastic often used in thermoformed food and pharmaceutical packaging containers.
Thermoforming comprises a sequence of interdependent operations and is characterized by being sensitive to the intrinsic properties of thermoplastics, namely the lower heat conduction and the deformation capability strongly dependent on temperature. In general, the thermoforming comprises: a heating stage, which aims to allow the sheet to acquire the required deformability; a sheet deformation stage in order to reproduce the contours of the piece and finally, a cooling stage, which allows the part to be extracted from the mold without distorting. In this way, the final performance of thermoformed products results from the sum of all actions that occur in these three main stages. Since there are processing variables associated with each of the three stages, including the material properties as a function of temperature, optimizing the thermoforming process is a complex task.
Generally, the optimization of thermoforming, like in the other real word optimization problems, consists of relating the effect of the operating variables of each stage with the to deal with the objectives to be optimized, as will be described in Section 3. Additionally, most of these works are based on experimental data, which makes the process unrealistic given the time required to do the experiments.
In addition to the general studies on the optimization of thermoforming, it is possible to find in the literature studies aiming at controlling the different stages of the process, namely the heating and forming stages.
Regarding the heating stage, the main objective is controlling the sheet temperature and the temperature fields developing during this step. This is very complex since the heating methods are indirect, that is, the sheet temperature is controlled by adjusting the heaters variables, and not directly. Furthermore, determining an optimum processing window in thermoforming process with the aim of achieving high quality parts is critical. In practice, the infrared heating stage is crucial since the final thickness distribution of the thermoformed part is closely related to the temperature of the sheet.
Wang and Nied [6] used a numerical approach based on the finite element method to obtain inverse solutions for the thermoforming processes. This was done by specifying a desired final thickness distribution and iteratively solving the system for the temperature field needed to obtain the desired result. A uniform initial temperature distribution is used as the initial guess for the iterative optimization procedure. Subsequently, updated non-uniform initial temperatures are obtained, based on the complete process of simulation for achieving the final thickness distribution during thermoforming. Suitable inverse solutions are achieved once the desired thickness distribution is obtained within a specified tolerance. The sensitivity of final part thickness to perturbations in temperature distribution is also investigated and shown to be a potential problem for precise thickness control in industrial applications.
Bordival et al. [7] developed an automatic optimization method of the ovens geometric parameters to be used in thermoforming. The first time, a simple analytical model, coupled to a nonlinear constraint optimization method (Sequential Quadratic Programming) allows one to find the best set of parameters, according to a cost function representing, for example, the heat flux uniformity. Then, with these optimized parameters, an accurate raytracing method is used to compute the irradiation resulting from the interaction between lamps and the thermoplastic sheet. Finally, a control volume method is implemented to solve the threedimensional transient heat transfer equation, where the radiation source is approximated by a diffusion Rosseland model.
Chy and co-authors [8,9] presented a method to control the surface temperature of a plastic sheet using model predictive control (MPC) to solve the inverse heating problem (IHP). The model was implemented on a complex thermoforming oven with a large number of inputs and outputs for precise control of sheet temperatures under hard constraints on heater temperature and their rates. Even though the MPC controller can handle a multivariable process, the large number of computations makes it difficult to apply to large systems such as multi-zone temperature control in a thermoforming machine.
Li and co-authors [10,11] suggested a methodology to compute and optimize the sheet temperature by controlling the heater temperature. The steady-state optimum distribution of heater power is first ascertained by a numerical optimization to obtain a uniform sheet temperature. The time-dependent optimal heater input is then determined to decrease the temperature difference through the direction of the thickness using the response surface method and the D-optimal method. The results show that the time-dependent optimum heater power distribution gives an acceptable uniform sheet temperature in the forming temperature range by the end of the heating process.
More recently, Erchiqui and co-authors [12][13][14][15] proposed the application of two different meta-heuristic algorithms, Simulated Annealing (SA) and Evolutionary Algorithms (EA) to detect, from a fixed and random set of temperatures of the radiant zones of the oven, the best temperatures that must be assigned to the heating zones in order to ensure a uniform sheet temperature. For numerical heating analysis, the nonlinear heat conduction problem is solved by a specific 3D volumetric enthalpy-based computational method.
However, these studies only considered the effects of the heating phase on the final thickness distribution of the thermoformed parts, having a starting point a sheet with constant thickness, that is, the aim is to determine how to heat the different zones of the sheet in order to induce uniform thickness on the final part.
Furthermore, an important limitation of the methodologies proposed in the literature lies in the fact that the process is intrinsically multi-objective, and the different stages cannot be considered independent. In reality, the results of one stage of the process depend strongly on the remaining stages, the final part thickness distribution does not depend only on the heating, but also on the initial thickness of the sheet used. Therefore, the main aim of the work is to apply multi-objective strategies to optimize the plastics thermoforming process. The intention is to propose a more global methodology that can take into account the different steps and characteristics of this process, as described in the next section. However, due to its complexity, only the stage of the sheet deformation will be considered here. More specifically, the aim is to study the influence of the initial sheet thickness distribution on the optimization of the process. For that purpose, three different initial thickness distribution sheets will be considered. As far as we know, there is no work in the literature that addresses the problem from this point of view. From an industrial perspective, the fabrication of the sheets to be thermoformed must be changed if economic and/or environmental gains can be obtained.
The paper is organized as follows: in Section 2 the details of thermoforming related to optimization are explained, Section 3 addresses the concepts of multi-objective optimization and the algorithm used, in Section 4 the results and discussion for the optimization of thermoforming are presented and in Section 5 the conclusions are stated. Figure 1 illustrates schematically the thermoforming phases and the optimization sequence. The phase order is indicated by open arrows and consist basically of: (A) producing a plastic sheet, typically by extrusion; (B) heating the sheet until it can be deformed without breaking; (C) shaping the sheet against the contours of a mold by applying a pressure difference on both sides, either by vacuum or pressurized air; (D) cooling the product obtained to make it possible to remove it from the mold and (E) remove the part from the mold. Each one of these phases has particular characteristics that must be considered when optimizing the process, either in terms of what concerns the operating conditions (e.g., heating and cooling times, air pressure and oven temperatures) and/or design parameters (e.g., sheet thicknesses, heater location, heating methods and mold geometry) [16,17]. Figure 2 shows schematically how the shaping process evolves. The heated sheet is forced to deform by the action of a pressure differential between both sides of the sheet. Initially, this deformation is uniformly distributed, but when it touches the cold mold surfaces the plastic material is frozen and only the remaining part of the sheet continues to deform. This implies that the last part of the sheet touching the mold will present the lowest thickness, since that the total volume of the sheet is conserved. This corresponds to the region of the lower corners in Figure 2B.  Figure 2 shows schematically how the shaping process evolves. The heated sheet is forced to deform by the action of a pressure differential between both sides of the sheet. Initially, this deformation is uniformly distributed, but when it touches the cold mold surfaces the plastic material is frozen and only the remaining part of the sheet continues to deform. This implies that the last part of the sheet touching the mold will present the lowest thickness, since that the total volume of the sheet is conserved. This corresponds to the region of the lower corners in Figure 2B. Therefore, the most important objective when producing a thermoformed part is to guarantee that its final thickness is as uniform as possible. This requirement is necessary in order to be possible to accomplish two main objectives, mainly in the shaping phase: to minimize the amount of material necessary and to maximize the mechanical behavior and/or other required characteristics.    Figure 2 shows schematically how the shaping process evolves. The heated sheet is forced to deform by the action of a pressure differential between both sides of the sheet. Initially, this deformation is uniformly distributed, but when it touches the cold mold surfaces the plastic material is frozen and only the remaining part of the sheet continues to deform. This implies that the last part of the sheet touching the mold will present the lowest thickness, since that the total volume of the sheet is conserved. This corresponds to the region of the lower corners in Figure 2B. Therefore, the most important objective when producing a thermoformed part is to guarantee that its final thickness is as uniform as possible. This requirement is necessary in order to be possible to accomplish two main objectives, mainly in the shaping phase: to minimize the amount of material necessary and to maximize the mechanical behavior and/or other required characteristics.  Therefore, the most important objective when producing a thermoformed part is to guarantee that its final thickness is as uniform as possible. This requirement is necessary in order to be possible to accomplish two main objectives, mainly in the shaping phase: to minimize the amount of material necessary and to maximize the mechanical behavior and/or other required characteristics.

The Process
To produce a part with uniform thickness, it is possible to act in different phases of the process. First, as shown in Figure 3, it is possible to produce sheets with different thicknesses at different regions of the part to be produced. Therefore, in the regions where more deformation occurs during the shaping phase, and, as a consequence, the thickness of the part will be smaller, it is possible to increase the thickness of the original sheet. Three situations are possible, as illustrated in Figure 3: (A) constant sheet thickness, (B) variable sheet thickness in one direction (x) and (C) variable sheet thickness in two directions (x and y). However, the processes used to produce the sheets are different, which must be taken into account in any optimization process since the costs can be very different as well. While (A) and (B) can be produced by extrusion, (C) must be produced by injection molding, given the specific and complicated geometry. While the sheet with constant thickness (A) can be produced in the usual way, using a flat die and a calender, to produce sheet (B) it is necessary to change the die and the calender rolls in order to induce different thicknesses.
Three situations are possible, as illustrated in Figure 3: (A) constant sheet thickness, (B) variable sheet thickness in one direction (x) and (C) variable sheet thickness in two directions (x and y). However, the processes used to produce the sheets are different, which must be taken into account in any optimization process since the costs can be very different as well. While (A) and (B) can be produced by extrusion, (C) must be produced by injection molding, given the specific and complicated geometry. While the sheet with constant thickness (A) can be produced in the usual way, using a flat die and a calender, to produce sheet (B) it is necessary to change the die and the calender rolls in order to induce different thicknesses. Secondly, the sheet deformation during the shaping phase depends on the mechanical properties of the polymeric material used at the heating temperature. This temperature dependency makes it possible to obtain different sheet deformations in order to control the thickness distribution, that is, in the regions where more deformation is required, the temperature must be lower [11,18]. Figure 4 illustrates schematically the process of heating the plastic sheet. Several parameters can be considered: (i) the number of heaters, for example, for sheets with higher thickness two heaters can be used, one each side; (ii) the distance between the heaters bank and the sheet; (iii) the dimensions of each Secondly, the sheet deformation during the shaping phase depends on the mechanical properties of the polymeric material used at the heating temperature. This temperature dependency makes it possible to obtain different sheet deformations in order to control the thickness distribution, that is, in the regions where more deformation is required, the temperature must be lower [11,18]. Figure 4 illustrates schematically the process of heating the plastic sheet. Several parameters can be considered: (i) the number of heaters, for example, for sheets with higher thickness two heaters can be used, one each side; (ii) the distance between the heaters bank and the sheet; (iii) the dimensions of each individual heater and (iv) temperature of each individual heater. This approach is illustrated in Figure 4, where a mesh of heaters (e.g., ceramic heaters) can be used, each one with the possibility of having different temperature distributions. In any case, it is fundamental to control the correct temperature in the entire sheet surface. In practice, and given the lower heat conductivity of plastics, a temperature gradient along the thickness of the sheet will arise. Furthermore, and if radiation heating is used, due to geometric reasons, the temperature at the center of the sheet will be higher than in its edge [19,20]. From the practical point of view and after the heating phase, the temperature of the entire sheet must be in the range of the forming temperatures. This temperature window is intrinsic to each polymer.
sheet must be in the range of the forming temperatures. This temperature window is intrinsic to each polymer.
Third, a straightforward way of minimizing the difference of thicknesses in the part thermoformed is related to the process of shaping. For that purpose, several thermoforming variants can be used: female mold (as illustrated in Figure 2), male mold, plug assisted, bubble forming, bubble forming and plug assisted and other combinations of the above [21][22][23][24][25].

Numerical Modeling
During the inflation process, the polymer sheet deforms due to the pressure difference imposed on both sides. Additionally, due to the deformation process, internal stresses are generated, which will limit the deformation promoted by the pressure difference. If the balance from those two loads is not null the velocity of the sheet changes, that is, there is a non-null acceleration.
For the adopted numerical approach, the polymer sheet is discretized into different computational cells as illustrated in Figure 5A. Moreover, the sheet is assumed to be thin enough to allow resorting to a membrane formulation, thus it is discretized into planar triangular cells each one with thickness ( Figure 5B). Considering the conditions above, the conservation of linear momentum is given by: where is time, ρ is the density, is the displacement vector, is the Cauchy stress tensor and is the gravitational acceleration vector. Following the classical Finite Volume Method approach, Equation (1) is integrated into the space in a computational cell with volume Ω and surrounded by the surface Γ, which, after applying the Gauss Divergence theorem results in: Third, a straightforward way of minimizing the difference of thicknesses in the part thermoformed is related to the process of shaping. For that purpose, several thermoforming variants can be used: female mold (as illustrated in Figure 2), male mold, plug assisted, bubble forming, bubble forming and plug assisted and other combinations of the above [21][22][23][24][25].

Numerical Modeling
During the inflation process, the polymer sheet deforms due to the pressure difference imposed on both sides. Additionally, due to the deformation process, internal stresses are generated, which will limit the deformation promoted by the pressure difference. If the balance from those two loads is not null the velocity of the sheet changes, that is, there is a non-null acceleration.
For the adopted numerical approach, the polymer sheet is discretized into different computational cells as illustrated in Figure 5A. Moreover, the sheet is assumed to be thin enough to allow resorting to a membrane formulation, thus it is discretized into planar triangular cells each one with thickness tk ( Figure 5B).
Considering the conditions above, the conservation of linear momentum is given by: where t is time, ρ is the density, u is the displacement vector, σ is the Cauchy stress tensor and g is the gravitational acceleration vector. Following the classical Finite Volume Method approach, Equation (1) is integrated into the space in a computational cell with volume Ω and surrounded by the surface Γ, which, after applying the Gauss Divergence theorem results in: which shows that the acceleration of the computational cell is affected by the surface and body forces. For the body forces, only the weight due to gravity is considered, which can be obtained by: where V c is the computational cell volume. which shows that the acceleration of the computational cell is affected by the surface and body forces. For the body forces, only the weight due to gravity is considered, which can be obtained by: where is the computational cell volume. The surface forces are of two types, namely external and internal forces. Following the membrane formulation, the external forces are due to the pressure difference between the inner and outer surfaces and the internal forces result from the stress induced by the material deformation. Based on Equation (2), where, in each time step, the acceleration of each computational cell is an explicit function of the applied loads, the surface forces are the sum of the internal and external forces and can be calculated by the following equation: where is the imposed pressure difference, the internal stress vector, A the cell surface area (see Figure 5), with normal vector and A the cross-sectional cell area, with normal vector . The internal forces are calculated at each cell cross section along the thickness. Due to the quick deformation process that takes place in thermoforming, which has been shown to be mostly elastic since the original shape of the sheet if the plastic part is heated [26], for this problem, the relation between the internal stresses was assumed to be based on the elasticity theory: where ℂ is the stiffness fourth-order constitutive tensor and the strain tensor. In this, and must be work conjugate pairs of stresses and strains respectively. Considering the membrane formulation, the elastic solid undergoes through large displacements, large rotations and large strains and non-linearity cannot be neglected, the proper work conjugate pair of stresses and strains are the Second Piola-Kirchhoff stress tensor, for and the Green-Lagrange strain tensor, for , which is given by: The surface forces are of two types, namely external and internal forces. Following the membrane formulation, the external forces are due to the pressure difference between the inner and outer surfaces and the internal forces result from the stress induced by the material deformation. Based on Equation (2), where, in each time step, the acceleration of each computational cell is an explicit function of the applied loads, the surface forces are the sum of the internal and external forces and can be calculated by the following equation: where p is the imposed pressure difference, τ the internal stress vector, A S the cell surface area (see Figure 5), with normal vector n s and A f the cross-sectional cell area, with normal vector n f . The internal forces are calculated at each cell cross section along the thickness. Due to the quick deformation process that takes place in thermoforming, which has been shown to be mostly elastic since the original shape of the sheet if the plastic part is heated [26], for this problem, the relation between the internal stresses was assumed to be based on the elasticity theory: τ = C : where C is the stiffness fourth-order constitutive tensor and the strain tensor. In this, τ and must be work conjugate pairs of stresses and strains respectively. Considering the membrane formulation, the elastic solid undergoes through large displacements, large rotations and large strains and non-linearity cannot be neglected, the proper work conjugate pair of stresses and strains are the Second Piola-Kirchhoff stress tensor, for τ and the Green-Lagrange strain tensor, for , which is given by: In Equation (6), F is the deformation gradient tensor, given by the derivative of each component of the deformed position x vector with respect to each component of the reference position X vector, and I is the unit tensor. This formulation is equivalent to the hyperelastic Saint Venant-Kirchhoff model [27].
On the implemented explicit numerical procedure, at each time step all the loads are calculated with Equations (3)-(7) and the acceleration of each cell is given by Equation (2). This allows one to calculate the new cell position and advance to the next time step.
On all the numerical studies performed in this work an elastic modulus of 10 MPa was employed, which is representative of HIPS at the usually employed processing temperature of 140 • C [28]. Moreover, the material was assumed to be incompressible, thus, a Poisson's ratio of 0.5 was considered.

Experimental Assessment
The modeling computational results were assessed using the circular cup illustrated in Figure 6 [29]. The experimental work was performed on a modified Illig thermoforming machine that allows individual monitoring and control of the main thermoforming variables. The heater bank included 110 Elstein 125 W ceramic radiators (each 120 × 60 mm). Square portions of the sheet were cut and manually clamped (fixed) on all four sides, in such a way that the material was not allowed to slide into the mold cavity. Once the heating stage is completed and a uniform sheet temperature of 150 • C is reached, the frame (holding the sheet) moves toward the forming station. The lower mold section moves vertically, clamping the sheet circularly. The parts were vacuum-formed at 0.7 bar. After forming and cooling, the cups were extracted, cut into four equal portions along a line passing by the center of the base midpoint and their thicknesses along that line and circumferentially were measured with a caliper (precision of 0.01 mm). At least five cups were produced for each set of variables. The results are generally presented along the arc length but correspond to averages in the circumferential direction. The material used in the sheets is a mixture of 50% Polystyrene (PS Laquerene 1540) with 50% Impact Polystyrene (HIPS Laquerene 7240), both produced by Elf-Atochem.

Multi-Objective Evolutionary Algorithms
In a Multi-Objective Optimization Problem (MOP), the goal is to minimi objectives simultaneously, that is, to find feasible solutions where every objective fun is minimized. A multi-objective optimization problem with objectives and de variables can be formulated as follows: The theoretical results were computed using a commercial software T-SIM using a K-BKZ viscoelastic model whose material properties are described in [29] and by using the OpenFOAM considering the hyperelastic Saint Venant-Kirchhoff model for the material. As can be seen, the computational results show a good agreement with the experimental measurements, the behavior of both theoretical models being very similar. In fact, regarding the simulation of the thickness distribution in thermoforming, it is possible to find, throughout the literature, several works that use hyperplastic models (e.g., Ghobadnam [30], Bernard et al. [31], Oueslati et al. [32], Jeet et al. [33]) and viscoelastic models (e.g., Cha et al. [34], Atami et al. [28]). Although researchers are increasingly turning to viscoelastic models for thermoforming applications, according to O'Connor et al. [35], there is little agreement on the best model to use for any given polymer material.
Furthermore, the use of hyperelastic models is based on the experience reported by Schmidt and Carley [26] and Rosenzweig et al. [36]. They described the "plastic memory", where blown bubbles returned completely to their original flat sheet form, either by suddenly releasing the forming pressure or by annealing spheroidal bubbles at a proper temperature. The observation of complete recovery has led the authors to conclude that elastic-like behavior can be ascribed to the processes of blowing soft plastic sheets.
Therefore, and considering that the main objective of this work is to evaluate the optimization process and not to predict more accurately the thickness distribution, and for computation time reasons, a hyperelastic model was chosen to predict the plastic sheet behavior during the forming process.

Multi-Objective Evolutionary Algorithms
In a Multi-Objective Optimization Problem (MOP), the goal is to minimize all objectives simultaneously, that is, to find feasible solutions where every objective function is minimized. A multi-objective optimization problem with m objectives and n decision variables can be formulated as follows: where x is the decision vector, i.e., x ∈ Ω ⊆ R n , f is the objective vector of m objective functions, that is, f ∈ Ω ⊆ R m , g i are p inequality constraint functions and h j are q equality constraint functions and l and u are the vectors of the lower and upper bounds on decision variables, respectively. Solutions are compared in terms of Pareto dominance, that is, a solution x is said to dominate a solution y (x ≺ y), if and only if f i (x) ≤ f i (y), for all i ∈ {1, . . . , m} and f j (x) < f j (y) for at least one j ∈ {1, . . . , m}. If all objectives do not conflict with each other, there exists a unique solution that minimizes all the objectives. In general, there are multiple conflicting objectives giving rise to a set of optimal solutions-the Pareto optimal set, instead of a single optimal solution. A feasible solution, x * is Pareto optimal if and only if there is no other solution y ∈ Ω, y = x * , that y ≺ x * . These solutions are incomparable to each other since none of these solutions can be said to be better than others.
In order to facilitate the decision-making process, an a posteriori method in which the search for an approximation (as close and diverse as possible) to the Pareto optimal set is performed before the decision-making process. Thus, the decision-maker can select the most suitable solution from this set according to their preferences. The Pareto optimal set provides valuable information with respect to the trade-offs between alternatives.
There is a large variety of approaches to solving MOPs. The designated classical (or traditional methods) are known as scalarization methods in which the MOP is reformulated and solved as a single objective optimization problem [37]. However, usually, scalarization functions involve several parameters that can be difficult to define in order to obtain different approximations to the Pareto optimal solutions. Other approaches are based on Evolutionary Algorithms (EAs). EAs are meta-heuristics that mimic the process of evolution of a population of individuals over generations. The fittest individuals to the environment will survive over time and therefore will be likely to reproduce. The offspring are generated by genetic operators, such as crossover and mutation and inherit the parent characteristics. Basically, in EAs, a population of individuals representing a candidate solution to the problem evolves using two mechanisms: Selection and Variation. The quality of each individual is measured using a fitness function that, in optimization problems, is related to the objective or objectives functions for single or multi-objective optimization problems, respectively. The selection mechanism guarantees that the best individuals have a higher probability of being selected for generating offspring. The variation mechanism provides the generation of new individuals by application of genetic operators.
Multi-Objective Evolutionary Algorithms (MOEAs) are EAs for solving MOPs. There are advantages to using this type of algorithms. They work with a population of candidate solutions, which makes it possible to approximate the entire Pareto optimal set in a single run. The performance of any MOEA is strongly related to the efficacy of its Selection mechanism that guides the search in the objective space, balancing convergence and diversity and also the variation mechanism that is responsible for the generation of offspring. A common approach to simulating natural selection in MOEAs consists of assigning fitness values to individuals in the population according to its quality for the MOP being solved. In terms of the type of fitness function used, MOEAs can be classified into three different types: dominance-, scalarizing-and indicator-based algorithms.
Dominance-based approaches calculate an individual's fitness on the basis of the Pareto dominance relation [38] or according to different criteria [39]. Scalarizing-based approaches [40] incorporate traditional mathematical techniques based on the aggregation of multiple objectives into a single parameterized function. Indicator-based approaches use performance indicators for fitness assignment; pairs of individuals are compared using some quality measure such as the hypervolume indicator [41,42]. The fitness value reflects the loss in quality if a given solution is removed [42,43].

SMS-EMOA
The multi-objective evolutionary optimization algorithm used in this work is based on the SMS-EMOA [32] and implemented in MATLAB. The outline of the SMS-EMO is given by Algorithm 1.

Algorithm 1 SMS-EMOA
% Initialize population at random with µ individuals k ← 0 Repeat q k+1 ← generate(P k ) % generate offspring by genetic operators P k+1 ← selection(P k ∪ {q k+1 }) % select µ best individuals k ← k + 1 Repeat until the stopping criterion is fulfilled The search starts from an initial population of µ individuals randomly generated satisfying the boundary constraints of the decision variables. In each generation, one parent is selected at random from the population and a single offspring is produced by application of a variation procedure. In this procedure, a Gaussian mutation with covariance matrix adaptation [44] is applied to the parent to produce a single offspring. Next, a deterministic selection procedure selects the µ best individuals to the next generation. The selection involves the non-dominated ranking of the population and the computation of the hypervolume contribution of each individual of the population.
The non-dominated sorting procedure is illustrated in Figure 7. First, the P k ∪ {q k+1 } individuals are ranked according to a non-dominated sorting procedure defining f fronts of sets of non-dominated individuals [28]. A rank is assigned to each front representing its level of domination. All solutions belonging to each front are incomparable. The first front F 1 contains the non-dominated solutions in P k ∪ {q k+1 }. The second front F 2 contains all non-dominated solutions in P k ∪ {q k+1 }\F 1 . This procedure is repeated until all solutions in P k ∪ {q k+1 } are included in a front. Any solution in front F i+1 is dominated by at least one solution of front F i for i ≥ 1.
individuals are ranked according to a non-dominated sorting procedure defining fronts of sets of non-dominated individuals [28]. A rank is assigned to each front representing its level of domination. All solutions belonging to each front are incomparable. The first front 1 contains the non-dominated solutions in ⋃{ +1 }. The second front 2 contains all non-dominated solutions in ⋃{ +1 }\ 1 . This procedure is repeated until all solutions in ⋃{ +1 } are included in a front. Any solution in front +1 is dominated by at least one solution of front for ≥ 1. Afterward, the hypervolume contribution [39] of each individual in the last front F = P k ⋃{q k+1 } is computed. Hypervolume definition guarantees that any non-dominated solution will not be replaced by a dominated solution since non-dominated solutions will have a higher hypervolume contribution than dominated ones. Hypervolume allows one to obtain a well-distributed set of solutions in the objective space as well as to guide the search toward the Pareto optimal front.
The hypervolume contribution computation is straightforward for problems with two objectives [41]. The approximations to the ideal vector ( * ) and the nadir vector ( ) computed in the current population are used to normalize objectives functions to the same order of magnitude in the interval [0,1]. The normalized objective function for the -th objective function is computed by: The solutions of the worst-ranked non-dominated front are sorted in ascending order according to the values of the first normalized objective function. A sequence that is additionally sorted in descending order is obtained since the solutions are mutually nondominated. Then, the hypervolume contributions of the solutions can be obtained by computing the rectangle area as illustrated in Figure 8. Given a sorted front with solutions = { 1 , 2 … , } , the hypervolume contribution of solution ( ( ) ) is computed by [39]: where = 2, … , − 1.
In Figure 8, it can be seen that the hypervolume contribution of solution b is inferior to c, as the area covered by b is smaller, which implies that this solution will not be selected in the next generation. Afterward, the hypervolume contribution [39] of each individual in the last front P k ∪ {q k+1 } is computed. Hypervolume definition guarantees that any non-dominated solution will not be replaced by a dominated solution since non-dominated solutions will have a higher hypervolume contribution than dominated ones. Hypervolume allows one to obtain a well-distributed set of solutions in the objective space as well as to guide the search toward the Pareto optimal front.
The hypervolume contribution computation is straightforward for problems with two objectives [41]. The approximations to the ideal vector (z * i ) and the nadir vector (z nad i ) computed in the current population are used to normalize objectives functions to the same order of magnitude in the interval [0,1]. The normalized objective function f norm i for the i-th objective function is computed by: The solutions of the worst-ranked non-dominated front are sorted in ascending order according to the values of the first normalized objective function. A sequence that is additionally sorted in descending order is obtained since the solutions are mutually nondominated. Then, the hypervolume contributions of the solutions can be obtained by computing the rectangle area as illustrated in Figure 8. Given a sorted front with r solutions F = {s 1 , s 2 . . . , s r }, the hypervolume contribution of solution s i (I(s i )) is computed by [39]: where i = 2, . . . , r − 1.
Mathematics 2021, 9, x FOR PEER REVIEW 13 of 21 The best individuals in terms of the domination ranking are selected to be the progenitors of the next generation. On the last front, the individual with the worst hypervolume contribution is discarded. This process is repeated until the stopping criterion is fulfilled.

The Optimization Problem to Solve
In the present study, the cup illustrated in Figure 9 was thermoformed with constant In Figure 8, it can be seen that the hypervolume contribution of solution b is inferior to c, as the area covered by b is smaller, which implies that this solution will not be selected in the next generation.
The µ best individuals in terms of the domination ranking are selected to be the progenitors of the next generation. On the last front, the individual with the worst hypervolume contribution is discarded. This process is repeated until the stopping criterion is fulfilled.

The Optimization Problem to Solve
In the present study, the cup illustrated in Figure 9 was thermoformed with constant temperature, a female mold and three types of sheets (as shown in Figure 3): constant thickness, linear spline variation and concentric spline variation. The aim was to determine the sheet thickness profile in order to: (i) minimize the initial sheet volume, as it implies less material use (f 1 ); (ii) minimize the minimum thickness found in the cells of the mesh used in the modeling calculations without hindering its mechanical behavior, as it is related with the capacity of the polymer sheet deformability, representing indirectly a measure of the thickness heterogeneity (f 2 ); and (iii) minimize the thickness heterogeneity, that is, the difference between the thickness of the part and a reference thickness, as defined by Equation (10) where M is the number of points located in a line defining the center of the cup in direction x, t 0 is a reference thickness defined by the user and t i are the thicknesses of the M points. This is a bound constrained multi-objective optimization problem with the following decision variables and objectives limitations (dimensions in meters): where x i is the sheet thickness of the constant thickness along the x direction. The thickness along the x direction is then imposed using a spline variation based on the 10 control points (see Figure 3), or the thickness of the five control points determining the concentric thickness variation, from the center to the border of the circle represented in Figure 3C. First, three bi-objective problems were considered (Cases 1 to 3), one for each case of sheet thickness variation and taking into account objectives f 1 and f 2 , with the aim of showing the effectiveness of the methodology in solving this problem. Then, two bi-objective problems using objectives f 1 and f 3 were considered (Cases 4 and 5, respectively with spline and concentric variations), the aim being to approach the solutions to a more realistic industrial situation. The multi-objective evolutionary algorithm used in this work is based on the SMS-EMOA. Considering the characteristics of the multi-objective optimization problem being solved, different configurations and search mechanisms can be adopted for SMS-EMOA. Thus, as described in Section 3.2, Gaussian mutation with covariance matrix adaptation was selected as the variation procedure, since the problem being solved comprises continuous variables. simultaneously, the hypervolume metric was chosen to provide well-distributed alternative solutions in the objective space. The configuration, including the parameter values of the algorithm, takes into account the considerable computational effort required to compute the objective function values. The population size was set to 20 individuals. The selection was carried out using a uniform distribution and variation was performed by the CMA evolution strategy operator [8], which is designed to work with real number representations. The maximum number of generations was set to 20.
continuous variables. simultaneously, the hypervolume metric was chosen to provide well-distributed alternative solutions in the objective space. The configuration, including the parameter values of the algorithm, takes into account the considerable computational effort required to compute the objective function values. The population size was set to 20 individuals. The selection was carried out using a uniform distribution and variation was performed by the CMA evolution strategy operator [8], which is designed to work with real number representations. The maximum number of generations was set to 20.

Results and Discussion
For Case 1, the sheet with constant thickness, both objectives (volume and thermoformed part minimum thickness) are in harmony, which means that the solutions of all generations are located in a line and the single optimal solution is the one more near the minimum of these two objectives. The same does not occur for the other two cases. Figure 10 shows the random initial population and the non-dominated solutions of the last generation for optimization Case 2, which comprises seven solutions. In this case, the sheet thickness is a spline generated from the 10 decision variables represented by black dots in the graphs, provided in Figure 11. It is clear that there is improvement along the 20 generations. The Pareto optimal front of the last generation presents some gaps between the solutions due to the thermoforming problem characteristics, as the location of the 10 points used to generate the sheet spline thickness are fixed and equally spaced along the x axis, which limits the search space.
The sheet and final part thickness profiles of solutions Ps1, Ps4 and Ps7 are illustrated in Figure 11. As can be seen in the graphs, the thickness profiles perpendicular to the spline ( Figure 11-left), when moving from solutions Ps1 to Ps7 the final part profile, are more uniform. From the point of view of the Decision-Maker, this seems to indicate that Ps7 is the best solution for practical purposes. This evidences the great advantage of performing this type of multi-objective optimization. Not only does the algorithm converge to better solutions, but it can also provide the Decision-Maker a set of solutions from which he can choose the best one to be applied in the real thermoforming practice, in this case, solution Ps7. Another characteristic of the solutions shown is that for x equal

Results and Discussion
For Case 1, the sheet with constant thickness, both objectives (volume and thermoformed part minimum thickness) are in harmony, which means that the solutions of all generations are located in a line and the single optimal solution is the one more near the minimum of these two objectives. The same does not occur for the other two cases. Figure 10 shows the random initial population and the non-dominated solutions of the last generation for optimization Case 2, which comprises seven solutions. In this case, the sheet thickness is a spline generated from the 10 decision variables represented by black dots in the graphs, provided in Figure 11. It is clear that there is improvement along the 20 generations. The Pareto optimal front of the last generation presents some gaps between the solutions due to the thermoforming problem characteristics, as the location of the 10 points used to generate the sheet spline thickness are fixed and equally spaced along the x axis, which limits the search space.     Figure 12 shows the Pareto optimal solutions for all the cases studied. As can be seen, in Case 3 (concentric spline), the optimization converges to five non-dominated solutions, identified by black dots. The final part thickness profile of four solutions for this case, Pc1, Pc2, Pc4 and Pc5 are represented in Figure 13 (solution Pc3 was not represented due to its similarity with solution Pc2). The black dots identify the location of the points used to generate the symmetrical concentric spline represented by a dashed line, the decision variables. Again, from solution Pc1 to solution Pc5, the profile obtained is more uniform, as in Case 2. Also, it is important to note that, in this case, the part thickness profile is the same in all directions, as the sheet thickness presents an axisymmetric distribution. The sheet and final part thickness profiles of solutions Ps1, Ps4 and Ps7 are illustrated in Figure 11. As can be seen in the graphs, the thickness profiles perpendicular to the spline (Figure 11-left), when moving from solutions Ps1 to Ps7 the final part profile, are more uniform. From the point of view of the Decision-Maker, this seems to indicate that Ps7 is the best solution for practical purposes. This evidences the great advantage of performing this type of multi-objective optimization. Not only does the algorithm converge to better solutions, but it can also provide the Decision-Maker a set of solutions from which he can choose the best one to be applied in the real thermoforming practice, in this case, solution Ps7. Another characteristic of the solutions shown is that for x equal to zero ( Figure 11-left) the sheet spline thickness is the same in all cases (t ≈ 2.0 × 10 −3 m), which produces identical part profiles in the transversal direction, as illustrated by the graphs of Figure 11-right. Figure 12 shows the Pareto optimal solutions for all the cases studied. As can be seen, in Case 3 (concentric spline), the optimization converges to five non-dominated solutions, identified by black dots. The final part thickness profile of four solutions for this case, Pc1, Pc2, Pc4 and Pc5 are represented in Figure 13 (solution Pc3 was not represented due to its similarity with solution Pc2). The black dots identify the location of the points used to generate the symmetrical concentric spline represented by a dashed line, the decision variables. Again, from solution Pc1 to solution Pc5, the profile obtained is more uniform, as in Case 2. Also, it is important to note that, in this case, the part thickness profile is the same in all directions, as the sheet thickness presents an axisymmetric distribution.    Figure 14 shows the part thickness profile of the unique solution found for Case 1, Pf1. As can be seen, the profile obtained is very similar to that of solutions Ps7 and Pc5 of the previous cases studied. This confirms that the optimization strategy is working and identify solutions with physical meaning, since different starting points, corresponding to diverse initial sheet thickness profiles, led to identical solutions when those solutions are located in the same region of the search space, as can be seen in Figure 12.   Figure 14 shows the part thickness profile of the unique solution found for Case 1, Pf1. As can be seen, the profile obtained is very similar to that of solutions Ps7 and Pc5 of the previous cases studied. This confirms that the optimization strategy is working and identify solutions with physical meaning, since different starting points, corresponding to diverse initial sheet thickness profiles, led to identical solutions when those solutions are located in the same region of the search space, as can be seen in Figure 12. Finally, Figure 15 shows the non-dominated solutions of the 20th generation for Cases 4 and 5, considering spline and concentric sheet thickness variation and t0 equal to 0.5 mm. As previously, the gaps between the solutions are due to the problem constraints related to a limited search space. It is clear that, as expected, the concentric variation produces much better results concerning the uniformity of the thickness. Solutions P's3 and P'c3 are the same as those obtained previously, that is, solutions Ps7 and Pc5. However, the solutions found depend strongly on the value chosen for t0. A practical alternative consists of finding the desired thickness profile for the final part through a mechanical analysis, for example, as suggested in the scheme of Figure 1.
Also, it is important to note that technically and economically it is difficult to produce concentric initial sheet variations, as it implies high costs, mainly when the parts are to be produced in big quantities, as is commonly the case.  Finally, Figure 15 shows the non-dominated solutions of the 20th generation for Cases 4 and 5, considering spline and concentric sheet thickness variation and t 0 equal to 0.5 mm. As previously, the gaps between the solutions are due to the problem constraints related to a limited search space. It is clear that, as expected, the concentric variation produces much better results concerning the uniformity of the thickness. Solutions P's3 and P'c3 are the same as those obtained previously, that is, solutions Ps7 and Pc5. However, the solutions found depend strongly on the value chosen for t 0 . A practical alternative consists of finding the desired thickness profile for the final part through a mechanical analysis, for example, as suggested in the scheme of Figure 1.  Also, it is important to note that technically and economically it is difficult to produce concentric initial sheet variations, as it implies high costs, mainly when the parts are to be produced in big quantities, as is commonly the case.

Conclusions
Faced with the complexity of the real optimization problem in the field of plastics thermoforming described in this paper, the applicability of a multi-objective optimization strategy was proposed to deal with the process forming phase. The aim was to determine the better sheet thickness distribution that allows the production of parts with the least amount of material while assuring the appropriate characteristics of the final part. A reduction of circa 30% in the volume of the material used can be obtained when comparing the solutions of the initial population and the ones in the Pareto front.
The main scientific contribution of this work is to add a new approach to control the thickness distribution of thermoformed parts, in addition to the variations of technique and differential heating previously mentioned. Since the final thickness of the parts depends on the initial thickness of the sheet, several solutions are proposed for adjusting the initial thickness of the sheet in order to achieve the most favorable thickness.
The results obtained showed that the methodology proposed was able to capture the singular features of the process allowing us to conclude that the strategy proposed might be successfully applied in the optimization of the various steps of plastics thermoforming. This work constitutes a further step to support design approaches associated with this important plastics processing technology.