A metamodel-based multicriteria shape optimization process for an aerosol can

1. Introduction

In the wide range of metallic packaging, the aerosol can is one of the most commonly used. One of the main requirements of specifications asked by the customers is the resistance to internal pressure. In fact, the aerosol cans are filled with fluid at high pressure; for this reason, the structural stability of their bottom is then delicate to maintain. The goal of the canmakers is to look for new profiles of the aerosol can’s bottom that resists against the internal pressure of the fluid (Fig. 1).

In the present work, we address the problem of shape optimization of the can’s bottom, in order to control the dome growth DG (e.g., displacement of can base) at a proof pressure as well as the dome reversal pressure DRP, a critical pressure at which the aerosol can’s bottom looses stability (e.g., initiates buckling). Those two criteria are known to be conflicting; therefore, our aim was to identify the Pareto front of this problem [1].

In most cases, the identification of the Pareto front for industrial optimization problems (e.g., all the time are structural problems) is very costly in computing time because we need a large number of evaluations for the criteria to be optimized. To overcome this, the most researchers and industry begin to develop algorithms which consist of coupling the conventional methods to capture the pareto front with metamodels aimed at cheap costs evaluation. There are two ways to do this coupling: The first idea is to lead optimization with a dedicated algorithm and use an updated metamodel for a certain number of evaluations until finding the solutions (e.g., strong coupling). The second idea is to lead optimization with the metamodel and only do the exact calculations of the obtained solutions (e.g., weak coupling). There are several descriptions of the coupling techniques, see for example [2], [3], [4], [5]. In this context, we can cite as examples, some research with different applications as the automotive structures [6], [7], [8], [9], aerodynamic design [10], [11], [12]and the building energy [13].

For our work, we have used one of the known methods adapted to the approximation of pareto front which is the NNCM Method [14], [15], [16]coupled a weak coupling to the RBF metamodel [17], [18], [19]. Firstly, we have studied the consistency of the NNCM-RBF coupling by solving several standard mathematical benchmarks; then, the coupling is applied to our industrial case, namely the shape optimization of the can’s bottom.

The structure of this paper is as follows: Section 2 summarizes the methods used to develop our algorithm with a validation for some academic test cases and Section 3 demonstrates the industrial case with the obtained results. Conclusions and future contributions are listed in Section 5.

2. NNCM coupled with the RBF metamodel

This section presents the coupling of NNCM method and the RBF metamodel and its validation through several test cases.

2.1. NNCM

The normalized normal constraint method is an approach proposed by messac and Yehya for generating a set of evenly spaced solutions on a Pareto frontier - for multicriteria optimization problems - [14].

Throughout our work, we address a special case of multicriteria optimization problems (two objective functions), and we define the mathematical representation of this optimization problem as follows:(1) $\begin{matrix} (P) \\ \{\begin{matrix} \min_{x} F (x) = (f_{1} (x), f_{2} (x)) \\ subject to \{x \in (C) \end{matrix} \end{matrix}$

In this formulation, the set $(C)$ contains all three types of constraints (equality, inequality and an upper and a lower bound constraints). The points x that fulfill all the constraints are feasible, while all other points are unfeasible.

For the NNCM method, the optimization takes place in the normalized objective space, and the main steps of the method can be summarized in the following points:

Let $x_{i}^{*}$ the respective global minimizers of $f_{i} (x)$ i = 1, 2, over $x \in (C)$ .

The Utopia and Nadir points are hypothetical points defined respectively by the best and the worst values for each of the criteria in the set (Fig. 2). They are written mathematically by the two following formulas $f^{U} = [f_{1} (x_{1}^{*}), f_{2} (x_{2}^{*})]$ and $f^{N} = [f_{1} (x_{2}^{*}), f_{2} (x_{1}^{*})]$ respectively.

The distances $L_{1}$ and $L_{2}$ , lying between the two global minimizer points and the Utopia point, are the elements of the following matrix:(2) $L = f^{N} - f^{U}$

Using the above definitions, the normalized objective space $\overline{F}$ can be evaluated as follows:(3) $\overline{F} (x) = ({\overline{f}}_{1} (x) = \frac{f_{1} (x) - f_{1} (x_{1}^{*})}{L_{1}}, {\overline{f}}_{2} (x) = \frac{f_{2} (x) - f_{2} (x_{2}^{*})}{L_{2}})$

Let us define ${\overline{N}}_{1}$ as the direction from $x_{1}^{*}$ to $x_{2}^{*}$ which is called Utopia line:(4) ${\overline{N}}_{1} = {\overline{F}}_{2}^{*} - {\overline{F}}_{1}^{*}$ where ${\overline{F}}_{i}^{*} = [{\overline{f}}_{1} (X_{i}^{*}), {\overline{f}}_{2} (X_{i}^{*})]$

We generate a set of evenly distributed points on the Utopia line as follows:(5) ${\overline{X}}_{pk} = \sum_{k = 1}^{2} β_{pk} {\overline{F}}_{k}^{*}$ where $\sum_{k = 1}^{2} β_{pk} = 1; 0 ⩽ β_{pk} ⩽ 1$ .

Using the set of evenly distributed points on the Utopia line, we generate a corresponding set of Pareto points by solving a succession of optimization runs of the $P_{NNCM}$ problem defined in (Eq. (6)). Each optimization run corresponds to a point on the Utopia line (Fig. 3).(6) $\begin{matrix} (P_{NNCM}) \\ \{\begin{matrix} \min_{x} {\overline{f}}_{2} (x) \\ subject to \{\begin{matrix} {\overline{N}}_{1} {(\overline{F} - {\overline{X}}_{pk})}^{T} ⩽ 0 \\ x \in (C) \end{matrix} \end{matrix} \end{matrix}$

2.2. RBF

Radial basis functions (RBF) have been developed for scattered multivariate data interpolation [19]. The method uses linear combinations of a radially symmetric function based on Euclidean distance or other such metric to approximate response functions. There exist many different kinds of radial basis function, and Schaback presents in his paper [20] a full Comparison of radial basis function interpolants. We conclude from this study that the simple classical form (Eq. (7)) remains the best and the most formula used by the majority of researchers:(7) $\tilde{f} (x) = \sum_{i = 0}^{n} ω_{i} ϕ (| | x - x_{i} | |)$ where n is the number of sampling points, x is the vector of design variables, $x_{i}$ is the vector of the $i$ th sampling point, $| | x - x_{i} | |$ is the Euclidean distance (e.g., The norm is usually Euclidean distance, although other distance functions are also possible.), $ϕ$ is a basis function (for example, Gaussian one $ϕ (r) = e^{- a_{f} r^{2}}$ where $a_{f}$ is the attenuation coefficient $(0 ≺ a_{f} ⩽ 1)$ ), and $ω_{i}$ is the unknown weighting coefficient which is obtained by solving the linear system:(8) $f = A . ω$ where $f = {[f (x_{1}), \dots, f (x_{n})]}^{T}$ and $A_{i, j} = ϕ (| | x_{i} - x_{j} | |)$ ( $i = 1, \dots, n; j = 1, \dots, n$ ).

The RBF metamodel formula (Eq. (7)) can be written as follows:(9) $\tilde{f} (x) = \sum_{i = 0}^{n} ω_{i} e^{- a_{f} | | x - x_{i} | |^{2}}$

To have a good metamodel at the level of precision, we must make a judicious choice for choosing the sampling points $x_{i}$ and the attenuation factor’s value $a_{f}$ . The choice of these elements influences directly the results.

2.3. Effect of the attenuation factor

The attenuation factor in radial basis function has a critical influence on the accuracy of the interpolation model. We illustrate this with a numerical example. The RBF metamodel is constructed for the quadratic-sine function (Eq. (10)) using the same sampling points data −10 equally spaced points in [0,2]- (Fig. 4).(10) $f (x) = x (1 - x) \sin (2 π x), x \in [0, 2]$

In his paper [23], Rippa discussed several empirical methods for choosing the best attenuation factor and he presents his technique called “rippa technique” to search the optimized value of this factor based on the leave-one-out validation [21].

Algorithm 1

leave-one-out validation - LOOV

1: Input: N the total number of RBF sampling points

x_{i}

and

a_{0}

the initial attenuation factor

2: For

n = 1 : N

x^{(n)} = \{x_{1}, x_{2}, \dots, x_{n - 1}, x_{n + 1}, \dots, x_{N}\}

▷

{\tilde{f}}^{(n)} (x, a) = \sum_{m = 1, m \neq n}^{N}

ω_{m}^{(n)} ϕ (| | x - x_{m} | |)

E_{n} = f (x_{n}) - {\tilde{f}}^{(n)} (x_{n}, a) = \frac{ω_{n}}{A_{n, n}^{- 1}}

▷

6: End For

E (a) = {[E_{1}, E_{2}, \dots, E_{N}]}^{T}

a_{f} = \underset{a}{argmin} | | E (a) | |

9: Output: The optimized attenuation factor

a_{f}

The computation of $E (a)$ requires the solution of N linear equations of order $(N - 1) * (N - 1)$ , and the total number of operations is of order $N^{4}$ which can be very expensive. An efficient algorithm is given by Rippa requires only one LU decomposition to compute the $E (a)$ , he demonstrated that:(11) $E_{n} = \frac{ω_{n}}{A_{n, n}^{- 1}}, n = 1 : N$ where $ω_{n}$ is the $n$ th element of the solution vector $ω$ and $A^{- 1}$ is the invertible matrix of $A$ defined in the (Eq. (8)).

2.4. Effect of the sampling points

To show this effect, we use the RBF metamodel to approximate the quadratic-sine function with different databases of sampling points (e.g. database generated by three different methods (Fig. 5) and from the obtained results which are presented in the figure (Fig. 6), we can conclude that the uniform distribution gives us the best approximation.

This example is one of the full comparative study with several function tests presented in my thesis [22]. This study led us to choose the uniform distribution as a method to generate the sampling points in order to construct the RBF metamodel.

2.5. The NNCM RBF coupling

Our aim was to couple the NNCM method and RBF metamodel in order to have a simple algorithm with a reasonable and reduced calculation time to solve multicriteria optimization problems. The NNCM-RBF coupling is presented in detail (see Algorithm 2) and tested for several optimization problems known as test problems, which are mathematical explicit functions (Hab3d, Fonseca and Tanaka problem).

Algorithm 2

The NNCM-RBF algorithm

1: Input: Let

(P)

the problem to be solved, N the total number of RBF sampling points and

a_{0}

the initial attenuation factor

\{\begin{matrix} \min_{x} F (x) = (f_{1} (x), f_{2} (x)) \\ subject to \{\begin{matrix} x \in (C) \end{matrix} \end{matrix}

x^{(N)} = \{x_{1}, x_{2}, \dots, x_{N}\}

▷ Uniform RBF sampling points

a_{f}

= RIPPA

(f_{1}, f_{2}, x^{(N)}, a_{0})

▷ The optimized attenuation factor by Rippa technique

({\tilde{f}}_{1}, {\tilde{f}}_{2})

= RBF

(f_{1}, f_{2}, x^{(N)}, a_{f})

▷ The approximate functions by RBF metamodel

x_{i}^{*} = \min_{x} {\tilde{f}}_{i} (x) i = 1 : 2

▷ Minimise each function

{\tilde{F}}^{u} = {[{\tilde{f}}_{1} (x_{1}^{*}), {\tilde{f}}_{2} (x_{2}^{*})]}^{T}

▷ Utopia point

{\tilde{F}}^{n} = {[{\tilde{f}}_{1} (x_{2}^{*}), {\tilde{f}}_{2} (x_{1}^{*})]}^{T}

▷ Nadir point

L = {\tilde{F}}^{u} - {\tilde{F}}^{n}

▷ Matrix of distences

{\overline{f}}_{i} (x) = \frac{{\tilde{f}}_{i} (x) - {\tilde{f}}_{i} (x_{i}^{*})}{L_{i}} i = 1 : 2

▷ Normalized functions

10:

{\overline{F}}_{i}^{*} = {[{\overline{f}}_{1} (x_{i}^{*}), {\overline{f}}_{2} (x_{i}^{*})]}^{T} i = 1 : 2

11:

{\overline{N}}_{1} = {\overline{F}}_{2}^{*} - {\overline{F}}_{1}^{*} i = 1 : 2

▷ Utopia line

12:

{\overline{X}}_{kj} = \sum_{k = 1}^{2} β_{k} {\overline{F}}_{k}^{*}, β_{i} ⩾ 0, \sum_{i = 1}^{2} β_{i} = 1

▷ generate points on the Utopia line

13: For each

β = {(β_{1}, β_{2})}^{T}

solve the

(P_{NNCM})

problem:

\{\begin{matrix} \min_{x} {\overline{f}}_{2} (x) \\ subject to \{\begin{matrix} {\overline{N}}_{1} {(\overline{F} - {\overline{X}}_{pk})}^{T} \\ x \in (C) \end{matrix} \end{matrix}

14: End For

15: Output: NNCM-RBF solutions

The results, (Fig. 7 and Table 1), show that the coupling NNCM-RBF converges to the Pareto frontier with fewer number (69%, 96% and 89%, for Fonseca, Tanaka and Hab3, respectively) of calls for the objective functions compared to a conventional NNCM.

Table 1. Functions call number required by NNCM and NNCM-RBF methods.

Problem	Method used	Prescribed Pareto points (P.N.)	Functions calls
Fonseca	NNCM	25	406
Fonseca	NNCM-RBF	25	128

Tanaka	NNCM	25	886
Tanaka	NNCM-RBF	25	32

Hab3	NNCM	25	704
Hab3	NNCM-RBF	25	64

After the validation of the coupling NNCM-RBF for the academic test cases, the practical case of shape optimization of the bottom of the aerosol can was considered.

3. Aerosol can

3.1. The aerosol can model

In the present study a 2D model provided by ArcelorMittal company and developed with ls dyna software was used to predict of the behavior of the aerosol can’s bottom when submitted to pressure. In fact, the aerosol can is considered axisymmetric; therefore, only a half generating two dimensional (2D) will be sufficient to represent the entire can. The assumption of 2D calculations is so near from the reality (Fig. 8).

For our work, we use two types of bottom’s cans (Fig. 9), the aerosol cans are made of thin high performance steel which have the characteristics shown in the Table 2.

Table 2. Process parameters used in simulation.

Process parameter

Value or law

Thick steel (e)

0.46 mm

Strain hardening exponent (n)

0.2

Yield strength (Re)

270 MPa

Ultimate strength (Rm)

380 MPa

Strength coefficient

K = e^{(Rm \cdot (n \cdot (1 - \ln (n))))}

Hollomon law

σ = K ∊^{n}

3.2. Calculation of criteria

During charging, there is pressure acting on the aerosol can’s bottom. To extract DRP value, it is necessary to calculate the pressure when the model deforms and it cannot resist more than threshold. This time is when the model switches from implicit to explicit calculation $T_{switch}$ . For extracting DG value, we need to just read the displacement of the point A on the foot of the aerosol can in the moment $T_{switch}$ (Fig. 10).

The figure (Fig. 11) shows an example on how to extract the criteria to be optimized in our industrial case by the ls dyna software.

We define the mathematical representation of the two criteria as follows:(12) $\{\begin{matrix} DRP = Pression (T_{switch}) \\ DG = Displacement ({(A)}_{T_{switch}}) \end{matrix}$

4. Shape optimization process

4.1. Shape optimization - problem description

The aim was to find an ideal new form of the bottom which satisfies the required bottom’s resistance according to the customer’s requirement better than the initial shape used which have the values given in the Table 3.

Table 3. Intial shapes of the aerosol can used - DRP and DG values.

Aerosol can	DRP (bar)	DG (mm)
N1	19.1092	0.8975
N2	15.2	0.4749

The goal is to figure out a design of the aerosol can’s bottom which satisfies a DRP value bigger than DRP of initial shapes and a DG value smaller than 1 mm.

Firstly, 4 points on the initial shape are selected as the design variables. By changing the position of each points (3 positions), a set of 81 points is obtained, where each point represents a given shape of the can bottom (Fig. 12). Then, for each point, the exact values of the two criteria DRP and DG are calculated. These values are collected to set a database of the sampling points allowing building the RBF metamodel for each criterion.

The design variables can only moving in one direction - in the vertical direction for $(φ_{1}, φ_{2}, φ_{3})$ and the horizontal one for $(φ_{1})$ .

Let $φ_{0} = (φ_{01}, φ_{02}, φ_{03}, φ_{04})$ the initial shape, and $α$ a positive offset. Then, we choose the bound constraints as follows:

$* φ^{lower} = (φ_{01} - α, φ_{02} - α, φ_{03} - α, φ_{04} - α)$ ,

$* φ^{upper} = (φ_{01} + α, φ_{02} + α, φ_{03} + α, φ_{04} + α)$ .

Finally, the industrial problem can be as the following formula:(13) $\begin{matrix} \max_{φ = (φ_{1}, φ_{2}, φ_{3}, φ_{4})} DRP (φ) / \min_{φ = (φ_{1}, φ_{2}, φ_{3}, φ_{4})} DG (φ) \\ subject to \{φ^{lower} ⩽ φ ⩽ φ^{upper} \end{matrix}$

4.2. Shape optimization - result and discussion

For $α = 0.5 mm$ , we computed an approximate Pareto front for the DG/DRP costs using our developed algorithm the NNCM-RBF coupling with different prescribed numbers of Pareto points (P.N. = 6, 12, 24). The figures (Figure 13, Figure 14) show the Pareto front obtained for the two study cases by our coupling and their exact evaluation. The Table 4 and the Table 5 show the overall time and total number of exact or surrogate evaluations used for the aerosol bottoms N1 and N2 respectively.

Table 4. Time required for the different functions call - Aerosol bottom N1 - (∗∗∗) = (∗) + (∗∗).

(P.N.)	Total ${time}^{* * *}$	Objective function		approximated function
(P.N.)	Total ${time}^{* * *}$	Call number	Time required∗∗	Call number	Time required∗
6	3 h 50 min 12 s	87	3 h 22 min 56 s	90,124	25 min 16 s
12	4 h 13 min 51 s	93	3 h 42 min 012 s	91,212	31 min 39 s
24	4 h 35 min 55 s	105	4 h 03 min 09 s	91,868	32 min 46 s

Table 5. Time required for the different functions call - Aerosol bottom N2 - (∗∗∗) = (∗) + (∗∗).

(P.N.)	Total ${time}^{* * *}$	Objective function		approximated function
(P.N.)	Total ${time}^{* * *}$	Call number	Time required∗∗	Call number	Time required∗
6	4 h 26 min 44 s	87	4 h 01 min 26 s	77,924	25 min 26 s
12	4 h 44 min 00 s	93	4 h 14 min 58 s	79,112	29 min 02 s
24	5 h 28 min 02 s	105	4 h 52 min 17 s	80,008	35 min 45 s

We can notice from the tables (Table 4, Table 5) that our approach has allowed us to save a remarkable computational time. For example, if we take the case (P.N. = 24) from the second case (N2) (Table 5), there are 105 calls of exact function evaluations and 80,008 for approximated function, respectively, which represent 0.14% and 99.86% of the total function calls used in our approach. But at the same time, we note that only this 0.14% of total calls take 89.10%of the total computing time required. This last remark justifies why we chose not to apply roughly the NNCM method with exact evaluations to solve this industrial case.

A simple comparison between the results obtained by our approach and the accurate evaluation of these solutions (Figure 13, Figure 14) allows us to assess that our results remain good ones notwithstanding the complexity of our cases study. We remarked also that all solutions are almost located at the boundary of the space formed by the elements of the RBF database (e.g. sampling points) (Fig. 15). Then we can conclude that the solutions obtained are likely NNCM solutions and our approach is able to solve the industrial problem with a reasonable computation time.