*Article* **Study of Parameters in the Genetic Algorithm for the Attack on Block Ciphers**

**Osmani Tito-Corrioso1,\* , Miguel Angel Borges-Trenard <sup>2</sup> , Mijail Borges-Quintana <sup>3</sup> , Omar Rojas <sup>4</sup> and Guillermo Sosa-Gómez 4,\***


**Abstract:** In recent years, the use of Genetic Algorithms (GAs) in symmetric cryptography, in particular in the cryptanalysis of block ciphers, has increased. In this work, the study of certain parameters that intervene in GAs was carried out, such as the time it takes to execute a certain number of iterations, so that a number of generations to be carried out in an available time can be estimated. Accordingly, the size of the set of individuals that constitute admissible solutions for GAs can be chosen. On the other hand, several fitness functions were introduced, and which ones led to better results was analyzed. The experiments were performed with the block ciphers AES(*t*), for *t* ∈ {3, 4, 7}.

**Keywords:** genetic algorithm; cryptanalysis; AES(*t*); optimization; heuristics

## **1. Introduction**

There are several methods and tools that are used as optimization methods and predictive tools. Several heuristic algorithms have been used in the context of cryptography; in [1], the Ant Colony Optimization (ACO) heuristic method was used, and a methodology with S-AES block encryption was tested, using two pairs of plain encrypted texts. In [2], a combination of GA and ACO methods was used for cryptanalysis of stream ciphers. In [3–5], the possibilities of combining and designing these analyzes using machine learning and deep learning tools were shown. In [6–8], the methods of the Artificial Neural Network (ANN), Support Vector Machine (SVM), and Gene-Expression Programming (GEP) were used as predictive tools in other contexts.

The Genetic Algorithm (GA) is an optimization method used in recent years in cryptography for various purposes, mainly to carry out attacks on various encryption types. Some of the research conducted in this direction is mentioned next. In [9], the authors presented a combination of the GA with particle swarm optimization (another heuristic method based on evolutionary techniques); they called their method genetic swarm optimization and applied it to attack the block cipher Data Encryption Standard (DES). Their experimental results showed that better results were obtained by applying their combined method than by using both methods separately. The proposal presented in [10] provided a preliminary exploration of GA's use over a Permutation Substitution Network (SPN) cipher. The purpose of the scan was to determine how to find weak keys. Both works [9,10] used a known plaintext attack, i.e., given a plaintext *T* and the corresponding ciphertext *C*, one is interested in finding the key *K*. In [10], the fitness function evaluates

**Citation:** Tito-Corrioso, O.; Borges-Trenard, M.A.; Borges-Quintana, M.; Rojas, O.; Sosa-Gómez, G. Study of Parameters in the Genetic Algorithm for the Attack on Block Ciphers. *Symmetry* **2021**, *13*, 806. https://doi.org/ 10.3390/sym13050806

Academic Editors: Juan Alberto Rodríguez Velázquez and Alejandro Estrada-Moreno

Received: 8 April 2021 Accepted: 27 April 2021 Published: 5 May 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

the bitwise difference (Hamming distance) between *C* and the ciphertext of *T*, using a candidate for the key, whereas, on the contrary, in [9] the Hamming distance between *T* and the decryption of the ciphertext of *C* is measured. In [11], a ciphertext-only attack on simplified DES was shown, obtaining better results than by brute force. The authors used a fitness function that combined the relative frequency of monograms, digrams, and trigrams (for a particular language). Since the key length was very small, they were able to use this kind of function. The approach in [12] was similar to [11]; it used essentially the same fitness function, but with different parameters. It was also more detailed regarding the experiments and compared them concerning brute force and random search. For more details on the area of cryptanalysis using GAs, see [13–15].

As in all evolutionary algorithms, it is always a difficulty in the GA that, as the number of individuals in the space of admissible solutions grows, in this case, the set of keys, it is necessary to perform a greater number of generations in order to obtain the best results. It is clear that the greater the number of generations, the more time the algorithm consumes, so it is important to be able to estimate the time that may be necessary to execute a certain number of desired generations. On the other hand, it is necessary to analyze fitness functions that allow obtaining better results with the fittest individuals obtained.

Symmetry is omnipresent in the universe; in particular, it is present in symmetric cryptography, where the secret key is known for both authorized parts in the communication channel essentially by symmetry. We worked with block ciphers, an important primitive of symmetric cryptographic, where the key space (the population of admissible solutions for the GA in this case) is exponentially big, making it impossible in many cases to fully move in that space.

In the present work, the ideas to divide the key space that were started in [16,17] were followed. Both methodologies for dividing the key space allow the GA search space to be reduced over a subset of individuals. For this case, we studied the behavior of time and the introduction of various fitness functions. The structure of the work is as follows. In Section 2, the general ideas of the GA and two methodologies for partitioning the key space are presented; in Section 3, several parameters of the cryptanalysis for block ciphers using the GA are studied; in Section 3.1, the time it takes to execute a certain number of iterations is analyzed, so that a number of generations to be carried out in an available time can be estimated; and in Section 3.2, other fitness functions are proposed. Finally, Section 4 gives the conclusions.

#### **2. Preliminaries**

#### *2.1. The Genetic Algorithm*

The GA is a heuristic optimization method. We assume that the reader knows the general ideas of how the GA works; see Algorithm 1. In this section, we briefly describe the GA scheme used in this work.

#### **Algorithm 1** Genetic algorithm.

**Input:** *m* (quantity of individuals in the population), *F* (fitness function), *g* (number of generations).

**Output:** the individual with the highest fitness function as the best solution.


The individuals from the populations are elements of the key space taken as binary blocks. For **Crossover**, the crossing by two points was used, and the crossover probability was fixed to 0.6. The **Mutate** operation consisted of interchanging the values of the bits of at most three random components of the binary block with a mutation rate of 0.2. The values of 0.6 and 0.2 were fixed for all experiments, and the study of the incidence of the variation of these values in the behavior of the GA was not addressed in this paper. An individual *x* is better adapted than another *y*, if it has greater fitness, i.e., if *F*(*x*) > *F*(*y*). Fitness functions are studied in more detail in Section 3.2. For the specification of the GA to block ciphers, see Section 3 of [16].

#### *2.2. Key Space Partition Methodologies*

The methodologies introduced in [16,17] allow GAs to work on a certain subset of the set of admissible solutions as if it were the complete set. The importance of this fact is that it reduces the size of the search space and gives the heuristic method a greater chance of success, assuming that the most suitable individuals are found in the selected subset. Let F *k*1 2 be the key space of length *k*<sup>1</sup> ∈ Z>0. It is known that F *k*1 2 has cardinality 2 *<sup>k</sup>*<sup>1</sup> , and therefore, there is a one-to-one correspondence between F *k*1 2 and the range <sup>h</sup> 0, 2*k*<sup>1</sup> <sup>−</sup> <sup>1</sup> i . If an integer *k*<sup>2</sup> is set, (1 < *k*<sup>2</sup> ≤ *k*1), then the key space can be represented by the numbers,

$$q 2^{k\_2} + r\_\* \tag{1}$$

where *q* ∈ h 0, 2*k*1−*k*<sup>2</sup> <sup>−</sup> <sup>1</sup> i and *r* ∈ h 0, 2*k*<sup>2</sup> <sup>−</sup> <sup>1</sup> i . In this way, the key space is divided into 2 *<sup>k</sup>*1−*k*<sup>2</sup> blocks (determined by the quotient in the division algorithm dividing by 2 *<sup>k</sup>*<sup>2</sup> ), and within each block, the corresponding key is determined by its position, which is given by the remainder *r*. The main idea is to stay in a block (given by *q*) and move within this block through the elements (given by *r*) using the GA. Note in this methodology that first *q* is set to choose a block, and then, *r* varies to be able to move through the elements of the block; however, the complete key in F *k*1 2 is obtained from Expression (1). We refer to this methodology as BBM. For more details on the connection with GAs, see [16].

The following methodology is based on the definition of the quotient group of the keys *G<sup>K</sup>* whose objective is to make a partition of F *k*1 2 in equivalence classes. It is known that F *k*1 2 , as an additive group, is isomorphic to Z<sup>2</sup> *k*1 . Let *h* be the homomorphism defined as follows:

$$\begin{array}{ccccc}\text{In}: \mathbb{Z}\_{2^{k\_1}} & \longrightarrow & \mathbb{Z}\_{2^{k\_2}}\\ & n & \longrightarrow & n \pmod{2^{k\_2}}\text{.}\end{array} \tag{2}$$

where *k*<sup>2</sup> ∈ Z><sup>0</sup> and 0 < *k*<sup>2</sup> < *k*1. We denote by *N* the kernel of *h*, i.e.,

$$N = \{ \mathbf{x} \in \mathbb{Z}\_{2^{k\_1}} | h(\mathbf{x}) = \mathbf{0} \in \mathbb{Z}\_{2^{k\_2}} \}. \tag{3}$$

Then, by the definition of *h*, we have that *N* is composed by the elements of Z<sup>2</sup> *k*1 , which are multiples of 2 *<sup>k</sup>*<sup>2</sup> . It is known that *N* is an invariant subgroup; therefore, the main objective is to calculate the quotient group of Z<sup>2</sup> *k*1 by *N*, and in this way, the key space will be divided into 2 *<sup>k</sup>*<sup>2</sup> equivalence classes. We denote by *G<sup>K</sup>* the quotient group of Z2 *k*1 by *N* (*G<sup>K</sup>* = Z<sup>2</sup> *k*1 /*N*). By Lagrange's theorem, we have that *o*(*GK*) = *o*(Z<sup>2</sup> *k*1 )/*o*(*N*), but *o*(*GK*) = *o*(Z<sup>2</sup> *k*2 ) = 2 *<sup>k</sup>*<sup>2</sup> , then,

$$o(N) = o(\mathbb{Z}\_{2^{k\_1}}) / o(\mathbb{Z}\_{2^{k\_2}}) = 2^{k\_1 - k\_2}.\tag{4}$$

Now, *N* can be described, taking into account that its elements are multiples of 2 *k*2 . For this, we take *<sup>Q</sup>* <sup>=</sup> {0, 1, 2, . . . , 2*k*1−*k*<sup>2</sup> <sup>−</sup> <sup>1</sup>}, then:

$$\begin{aligned} N &= \ &< \mathfrak{L}^{k\_2} > = \ & \{ \mathfrak{x} \in \mathbb{Z}\_{2^{k\_1}} \mid \exists \ \mathfrak{q} \in \mathbb{Q}, \ \mathfrak{x} = \mathfrak{q} \mathfrak{L}^{k\_2} \} \\ &= & \{ \mathfrak{0}, \mathfrak{L}^{k\_2}, \mathfrak{z} \ast \mathfrak{L}^{k\_2}, \mathfrak{z} \ast \mathfrak{L}^{k\_2}, \dots, (\mathfrak{L}^{k\_1 - k\_2} - 1) \ast \mathfrak{L}^{k\_2} \}. \end{aligned} \tag{5}$$

On the other hand,

$$G\_{\mathcal{K}} = \{N, 1+N, 2+N, \dots, (2^{k\_2}-2) + N, (2^{k\_2}-1) + N\}.\tag{6}$$

In this way, Z<sup>2</sup> *k*1 is divided into a partition of 2 *<sup>k</sup>*<sup>2</sup> classes given by *N*. *G<sup>K</sup>* is called the quotient group of keys. Let,

$$E: \{0, 1\}^m \times \{0, 1\}^n \to \{0, 1\}^n, m, n \in \mathbb{N}, m \ge n,\tag{7}$$

be a block cipher, *T* a plaintext, *K* a key, and *C* the corresponding ciphertext, i.e., *C* = *E*(*K*, *T*); *K* 0 is said to be a consistent key with *E*, *T*, and *C*, if *C* = *E*(*K* 0 , *T*) (see [16]). The idea here is also to go through, from the total space, the elements that are in a class and then find one (or several) consisting of the keys of that class. To be able to go through the elements of each class, note that Z<sup>2</sup> *k*2 is isomorphic with *GK*, and the isomorphism corresponds to each *r* ∈ Z<sup>2</sup> *k*2 its equivalence class *r* + *N* in *GK*; thus, selecting a class is setting an element *r* ∈ Z<sup>2</sup> *k*2 . On the other hand, the elements of *N* are of the form *q* 2 *k*2 (*q* ∈ *Q*); therefore, the elements of the class *r* + *N* are of the form,

$$q \, 2^{k\_2} + r \, \, q \in Q. \tag{8}$$

Then, the problem of looping through each element of each equivalence class consists of first setting an element of Z<sup>2</sup> *k*2 and then looping through each element of the set *Q*, to find a key of *G<sup>K</sup>* using Equation (8). The elements of the set *Q* have block length *k<sup>d</sup>* = *k*<sup>1</sup> − *k*2, and each class has 2 *<sup>k</sup><sup>d</sup>* elements. We refer to this methodology as TBB. Note that the TBB methodology is a kind of dual idea with respect to the BBM methodology, i.e., one first stays in the same class (given by *r*) and then moves within this class through the elements (given by *q*) using the GA. In this case, the length of the blocks is 2*k<sup>d</sup>* instead of 2*k*<sup>2</sup> .

The main difficulty in these methodologies is the choice of *k*2, since it is the parameter that determines the number of equivalence classes and, therefore, the number of elements within them. If in *GK*, *k*<sup>2</sup> increases, the classes have fewer elements, but there are more classes; on the contrary, if it decreases, so does the number of classes, but the number of elements of each increases. Something similar happens in the first methodology. The operations of the space partitioning and going through the elements of each class are done with the decimal representation and the specific operations of the GA with the binary representation. For more details, see [16,17].

In Figure 1, the relationship of the content by subsections and the attack on block ciphers are shown in a flowchart.

**Figure 1.** Flowchart of the relationship between content by subsections and the attack on block ciphers.

#### **3. Study of Parameters in the GA**

#### *3.1. Time Estimation*

In GAs, less complex operations such as mutation and crossing are performed within each class, where the elements have block length *k*<sup>2</sup> ≤ *k*<sup>1</sup> or *k<sup>d</sup>* ≤ *k*<sup>1</sup> depending on the way of partitioning the space. However, despite the variation of these two parameters, the calculation of the fitness function, being the function of greater complexity within the GA, is carried out using (8), i.e., with the complete key of length *k*1, and not with the part of it found in the class. This means that a variation in the number of elements in a class does not affect the fitness function's cost. Moreover, if all the parameters remain the same, the GA's time in each generation must be quite similar, even if *k*<sup>2</sup> varies. To check this, experiments were done with a PC with an Intel(R) Core (TM) i3-4160 CPU @ 3.60GHz (four CPUs), and 4GB of RAM. AES(*t*) encryption was used, a parametric version of AES, where *t* ∈ {3, 4, 5, 6, 7, 8} and also AES(8) = AES (see [18,19]). The experiment consisted of executing the GA with the BBM methodology and measuring the time (in minutes) that it took in a generation for different values of *k*<sup>2</sup> (keeping the other parameters fixed), then verifying if these data were used to forecast the time it would take in *n* generations. The size of the population was *m* = 100 in all cases.

Tables 1–3 summarize the results corresponding to AES(3), AES(4), and AES(7), respectively. The first column has the different values that were given to *k*2. The second column is the average time *tk*<sup>2</sup> that was obtained for a generation in 10 executions of each *k*2. The general mean for all the *k*<sup>2</sup> values is *t<sup>m</sup>* = 0.0435571 minutes approximately in Table 1, *t<sup>m</sup>* = 0.0519393 in Table 2, and *t<sup>m</sup>* = 0.1900297 in Table 3. The third column represents the number of generations (*ng*). The real-time that the algorithm takes, *t<sup>r</sup>* , appears in the fourth column. The fifth column is the estimated time, *t<sup>e</sup>* , that should be delayed, the calculation of which is based on:

$$t\_{\varepsilon} = t\_m n\_{\mathcal{S}}.\tag{9}$$

Finally, the last column is the error of the prediction, *E<sup>p</sup>* = |*t<sup>r</sup>* − *t<sup>e</sup>* |. With these experiments, we wanted to check for the procedure whether if for a specific value of *k*<sup>2</sup> and having *n<sup>g</sup>* generations, then the approximate time (*t*) that the GA would take to complete those generations was *t* ≈ *t<sup>e</sup>* .

With a generation, or very few, the average time it took for the GA was slightly slower, decreasing and tending to stabilize at a limit as it performed more iterations. This was due to probabilistic functions that intervened in the GA and a set of operations to randomly create an initial population. Therefore, the criterion for calculating the average time *tk*<sup>2</sup> was to let the GA finish executing in a certain number of generations, either because it found the key or because it reached the last iteration without finding it, and then calculate the average. Therefore, calculating *tk*<sup>2</sup> in a few generations or setting the amount to one, would get longer times; however, doing so would be valid if the intention were to go over the top in estimating the time that the algorithm consumed.


**Table 1.** Time estimation in AES(3).


**Table 2.** Time estimation in AES(4).

In the case of AES(7) (Table 3), we only experimented with the values 17 and 18 of *k*2, since considering all the previous (or higher) values would take a considerably longer time (given the greater strength of AES(7)).

**Table 3.** Time estimation in AES(7).


Similar results were obtained if more values of *k*<sup>2</sup> were chosen to calculate *tm*. For example, using a PC Laptop with a processor: Intel (R) Celeron (R) CPU N3050 @ 1.60GHz (two CPUs), ∼1.6 GHz, and 4 GB of RAM and going through all the values of *k*<sup>2</sup> from 10 to 48 (AES(3) key length), *t<sup>m</sup>* = 0.2340212 was obtained. Now, for *n<sup>g</sup>* = 215, we had *t<sup>r</sup>* = 48.14715 and *t<sup>e</sup>* = *tmn<sup>g</sup>* ≈ 50.3145. In another test: *n<sup>g</sup>* = 150, *t<sup>r</sup>* = 34.9565, then *t<sup>e</sup>* ≈ 35.1032. Note that the PC used in this case had different characteristics and less computational capacity than the experiments in Tables 1–3. The interesting thing is that under these conditions, the results were as expected as well.

In a similar way, the GA was executed with the TBB methodology for the search in *GK*, for values of *k<sup>d</sup>* equal to those of *k*<sup>2</sup> and different generations (*ng*). It was observed that the time estimates behaved in a similar way to the results presented previously for the BBM methodology. Note that in the AES(*t*) family of ciphers, the length of the key increases from 48 for AES(3) to 128 for AES(8); however, regardless of the key length, the same behavior was seen in all of them.

Now, we showed with these experiments another application of this study on time estimation. In the GA scheme with the BBM methodology, the total number of generations (iterations) to perform for a given value of *k*<sup>2</sup> is:

$$\mathbf{g} = \left\lfloor \frac{2^{k\_2}}{m} \right\rfloor. \tag{10}$$

Taking *n<sup>g</sup>* = *g*, by using *t<sup>e</sup>* , then we can do an a priori estimation for a given value of *k*2, of the total time it will take the GA to perform all the generations or a certain desired percent of them. For example, in AES(3), for *k*<sup>2</sup> = 16, in Expression (10), we have *g* = 655; now, since *t<sup>m</sup>* = 0.0435571 in Table 1, then the approximate time that the GA will consume to perform 655 generations is *t<sup>e</sup>* ≈ 655 · 0.0435571 ≈ 28.5299, as can be seen in the table. Another example can be seen in Table 2, also for *k*<sup>2</sup> = 16.

On the other hand, supposing we have an available time *t<sup>e</sup>* , to carry out the attack with this model, thus we may use (9) and (10), to compute an approximated value of *k*2, which implies doing the corresponding partition of the space and computing the number

of generations to perform for this time *t<sup>e</sup>* and the value of *k*2. In this sense, doing *n<sup>g</sup>* = *g* in (9), we have:

$$k\_2 \approx \left\lfloor \log\_2 \frac{mt\_\varepsilon}{t\_m} \right\rfloor. \tag{11}$$

We remark that the above is valid in the TBB methodology, only that *k<sup>d</sup>* is used instead of *k*2.

As can be observed, the results on the estimation of time were favorable. In this sense, the following points can be summarized:


#### *3.2. Proposal of Other Fitness Functions*

In the context of the BBM and TBB methodologies used in this work with the GA, we studied in this section which fitness functions provided a better response, in the sense that consistent keys were obtained as solutions in a greater percentage of occasions. Let *E* be a block cipher with length *n* of plaintext and ciphertext, defined as in Expression (7), *T* a plaintext, *K* a key, and *C* the corresponding ciphertext, that is *C* = *E*(*K*, *T*). Let:

$$D: \{0.1\}^m \times \{0.1\}^n \to \{0.1\}^n,\tag{12}$$

be the function of decryption of *E*, such that *T* = *D*(*K*, *C*). Then, the fitness function with which we have been working and based on the Hamming distance *dH*, for a certain individual *X* of the population, is:

$$F\_1(X) = \frac{n - d\_H(\mathbb{C}\_\prime E(X, T))}{n},\tag{13}$$

which measures the closeness between the encrypted texts *C* and the text obtained from encrypting *T* with the probable key *X* (see [16]). A similar function is the one that measures the closeness between plaintexts:

$$F\_2(X) = \frac{n - d\_H(T\_\prime D(X, \mathbb{C}))}{n}.\tag{14}$$

Another function that follows the idea of comparing texts in binary with *d<sup>H</sup>* is the weighting of *F*<sup>1</sup> and *F*2. Let *α*, *β* ∈ [0, 1] ⊂ R, such that *α* + *β* = 1, then this function would be defined as follows:

$$F\_3(X) = \mathfrak{a}F\_1(X) + \beta F\_2(X). \tag{15}$$

It is interesting to note that *F*<sup>3</sup> is more time consuming than each function separately, but the idea is to be more efficient in searching for the key.

The fitness functions proposed at this point are based on measuring the closeness of the plaintext and ciphertext, but in decimals. Let *Y<sup>d</sup>* be the corresponding conversion to decimals of the binary block *Y*. The first function is defined as follows,

$$F\_4(X) = \frac{2^n - 1 - |\mathcal{C}\_d - E(X, T)\_d|}{2^n - 1}. \tag{16}$$

Note that if the encrypted texts are equal, *C<sup>d</sup>* = *E*(*X*, *T*)*<sup>d</sup>* , then |*C<sup>d</sup>* − *E*(*X*, *T*)*<sup>d</sup>* | = 0, which implies that *F*4(*X*) = 1, i.e., if they are equal, then the fitness function takes the highest value. On the contrary, the greatest difference is the farthest they can be, i.e., *C<sup>d</sup>* = 2 *<sup>n</sup>* <sup>−</sup> <sup>1</sup> and *<sup>E</sup>*(*X*, *<sup>T</sup>*)*<sup>d</sup>* <sup>=</sup> 0, and therefore, *<sup>F</sup>*4(*X*) = 0. The following is a weighting of the functions *F*<sup>1</sup> and *F*4,

$$F\_5(X) = \alpha F\_1(X) + \beta F\_4(X). \tag{17}$$

Both functions have in common that they measure the closeness between ciphertexts. This is not ambiguous since, for example, if *C* and *E*(*X*, *T*) differ by two bits, the function *F*<sup>1</sup> will always have the same value no matter what these two bits are. On the contrary, it is not the same in *F*<sup>4</sup> if the bits are both more or less significant since the numbers are not the same in their decimal representation. The following function measures the closeness in decimals of plaintexts:

$$F\_6(X) = \frac{2^n - 1 - |T\_d - D(X, \mathbb{C})\_d|}{2^n - 1}. \tag{18}$$

Finally, the functions *F*7, *F*8, and *F*<sup>9</sup> are defined with respect to the previous ones as follows,

$$F\_7(X) \quad = \
u F\_2(X) + \beta F\_6(X),\tag{19}$$

$$F\_8(X) \quad = \ \ aF\_4(X) + \beta F\_6(X),\tag{20}$$

$$F\_{\overline{\mathbb{P}}}(X) \quad = \ a\_1 F\_1(X) + a\_2 F\_2(X) + a\_3 F\_4(X) + a\_4 F\_6(X),\tag{21}$$

where *α<sup>i</sup>* ∈ [0, 1] ⊂ R, *i* ∈ {1, 2, 3, 4} and 4 ∑ *i*=1 *α<sup>i</sup>* = 1. This guarantees that in general, each *Fj*(*X*) ∈ [0, 1] ⊂ R, *j* ∈ {1, 2, 3, 4, 5, 6, 7, 8, 9}.

The idea behind the introduction of these functions lies mainly in the fact that there are changes that the Hamming distance does not detect, as opposed to the decimal distance. For example, suppose the key is *a* = (1, 1, 1, 1, 1, 1)2, and *b* = (0, 0, 0, 0, 0, 1)<sup>2</sup> is the possible key, both in binary. It is clear that the Hamming distance is five, and the distance in decimals is 62 since *a* = 63 and *b* = 1; the fitness functions take the values 1 − 5/6 = 0.17 for the binary version and 1 − 62/63 = 0.016 for the decimal version. Now, if *b* = (0, 0, 1, 0, 0, 0)2, the binary fitness function would still be 0.17 since there are still five different bits; on the other hand, *b* = 8, so the decimal fitness function takes the value 1 − 55/63 = 0.13. Finally, if we take *b* = (1, 0, 0, 0, 0, 0)<sup>2</sup> = 32, then the distance in binary remains the same value, but the decimal continues to change, therefore, the fitness function as well, and takes the value 0.49. Therefore, this shows that the change of *b*, the decimal distance, is always detected, unlike the binary distance, which remains the same for certain changes.

AES(3) encryption attack experiments were carried out for the two methodologies for partitioning the key space to compare these functions. The main idea is to find the key and not do a component percent match analysis between them, where the fitness functions with the Hamming distance would be more useful. A PC with an Inter (R) Core (TM) i3-4160 CPU @ 3.60GHz (four CPUs), and 4 GB of RAM was used. For the results, we took into account the average time it took to find the key, the average number of generations in which it was found, the percentage of failures (in many attacks carried out), and a parameter called efficiency, *EF<sup>i</sup>* , which resulted in a weighting of the three previous criteria.

**Definition 1** (Fitness functions' efficiency)**.** *Let µ*1*, µ*2, *µ*<sup>3</sup> ∈ [0, 1] ⊂ R, *µ*<sup>1</sup> + *µ*<sup>2</sup> + *µ*<sup>3</sup> = 1*, tFi , i* = 1, · · · , *k, the time it takes the GA to find the key with F<sup>i</sup> , on an average for gF<sup>i</sup> generations,* *and pF<sup>i</sup> the percent of attempts in that the GA did not find the key with F<sup>i</sup> . Then, the efficiency, EF<sup>i</sup> , of the fitness function F<sup>i</sup> with respect to the other k* − 1 *functions, F<sup>j</sup>* , *j* 6= *i, is defined as,*

$$E\_{F\_i} = 1 - \left(\mu\_1 \frac{t\_{F\_i}}{\sum\_{\substack{k \\ \gamma = 1}}^k t\_{F\_\gamma}} + \mu\_2 \frac{g\_{F\_i}}{\sum\_{\substack{k \\ \gamma = 1}}^k g\_{F\_\gamma}} + \mu\_3 \frac{p\_{F\_i}}{\sum\_{\substack{k \\ \gamma = 1}}^k p\_{F\_\gamma}}\right). \tag{22}$$

Note that the number of generations and the failure percentage are inversely proportional to the efficiency *EF<sup>i</sup>* as the higher these parameters, the lower its efficiency fitness function. Table 4 presents the results of the comparison of the different fitness functions for the BBM space partitioning methodology, in this case *k* = 9. We took *α* = *β* = 0.5 and each *α<sup>i</sup>* = 0.25. To calculate *EF<sup>i</sup>* the values *µ*<sup>1</sup> = 0.33, *µ*<sup>2</sup> = 0.33 and *µ*<sup>3</sup> = 0.34 were taken for *tF<sup>i</sup>* , *gF<sup>i</sup>* , and *pF<sup>i</sup>* , respectively. Sorting *F<sup>i</sup>* with respect to efficiency, the first five would be *F*6, *F*8, *F*4, *F*5, and *F*2. It is noteworthy that of the first three that use only the Hamming distance, only *F*<sup>2</sup> appears.

**Table 4.** Comparison of fitness functions, with BBM.


In the comparison of these functions for the TBB methodology of partitioning the key space and searching in *GK*, the experiment results are presented in Table 5. In this case, ordering the functions by their efficiency, the first five would be *F*1, *F*4, *F*5, *F*8, and *F*6. Again, a single function appears from the first three, in this case *F*1, and the others repeat. Note in particular that *F*<sup>8</sup> (the weight of the functions in decimals) is better than *F*<sup>3</sup> (the weight of the functions in binary) in each of the parameters measured in both methodologies.

**Table 5.** Comparison of fitness functions, with TBB.


It is interesting to see what happens if the values of the weights are changed in the functions *F*5, *F*7, and *F*9, which combine the functions with distance in decimals and binary, keeping fixed *µ*1, *µ*2, and *µ*<sup>3</sup> for the calculation of *EF<sup>i</sup>* . In this sense, in the following group of experiments, the weights were assigned as follows for each methodology: the values

were 0.2 and 0.8; first, in each of these three functions, the subfunctions in binary were favored, from which *α* = 0.8, *β* = 0.2 (in *F*5, *F*7), *α*<sup>1</sup> = *α*<sup>2</sup> = 0.4, and *α*<sup>3</sup> = *α*<sup>4</sup> = 0.1 (in *F*9; note that this function has two subfunctions with the distance in binary and two in decimals); in this case, we identified the functions as *F*5*<sup>b</sup>* , *F*7*<sup>b</sup>* , and *F*9*<sup>b</sup>* ; then, we changed the order of these same weights, and the largest were given to the subfunctions whose distance was in decimals; and we identified the functions for this case as *F*5*<sup>d</sup>* , *F*7*<sup>d</sup>* , and *F*9*<sup>d</sup>* .

For the BBM methodology, the results are presented in Table 6. Note that according to *EFi* , the first is *F*7*<sup>d</sup>* , followed by *F*5*<sup>d</sup>* and *F*9*<sup>d</sup>* .


**Table 6.** Comparison of functions *F*5, *F*7, and *F*9, with BBM.

In Figure 2, these results are compared, according to *EF<sup>i</sup>* , with those of Table 4, also including the values of *F*5, *F*7, and *F*9. Sorting the functions according to their efficiency, the first five are *F*7*<sup>d</sup>* , *F*6, *F*8, *F*5*<sup>d</sup>* , and *F*9*<sup>d</sup>* .

**Figure 2.** Efficiency of all fitness functions in the BBM methodology.

Notice how the best results prevail in the functions with the distance in decimals. In this sense, *F*<sup>7</sup> and *F*<sup>9</sup> (now as *F*7*<sup>d</sup>* and *F*9*<sup>d</sup>* ) are incorporated into the first ones and three of those that already were in this group in the above experiments, *F*<sup>5</sup> (as *F*5*<sup>d</sup>* ), *F*6, and *F*<sup>8</sup>

In the case of the TBB methodology, the results are presented in Table 7. According to efficiency, the first is *F*7*<sup>b</sup>* , followed by *F*5*<sup>d</sup>* and *F*5*<sup>b</sup>* .

**Table 7.** Comparison of functions *F*5, *F*7, and *F*9, with TBB.


In Figure 3, these results are compared with those of all the functions of Table 5. The first five are now *F*1, *F*5, *F*4, *F*7*<sup>b</sup>* , and *F*5*<sup>d</sup>* ; notice how the functions that contain the distance prevail in decimals and this combined with binary. In the experiments, the best global behavior of the functions with the decimal distance is verified, and specifically in the BBM methodology, where the keys are grouped into intervals according to their decimal position in space, contrary to the other methodology, where the keys of each class are positioned throughout the space.

**Figure 3.** Efficiency of all fitness functions in the TBB methodology.

Note that when comparing Figures 2 and 3, the values of *EF<sup>i</sup>* that are in the tables are not directly compared, but rather, it is necessary to recalculate *EF<sup>i</sup>* taking into account that there are 15 functions. We mean,

$$E\_{F\_{\delta\_i}} = 1 - \left(\mu\_1 \frac{t\_{F\_{\delta\_i}}}{\sum\limits\_{\gamma=1}^k t\_{F\_{\delta\_\gamma}}} + \mu\_2 \frac{g\_{F\_{\delta\_i}}}{\sum\limits\_{\gamma=1}^k g\_{F\_{\delta\_\gamma}}} + \mu\_3 \frac{p\_{F\_{\delta\_i}}}{\sum\limits\_{\gamma=1}^k p\_{F\_{\delta\_\gamma}}}\right) . \tag{23}$$

where *δ<sup>i</sup>* ∈ {1, · · · , 9, 5*b*, 5*d*, 7*b*, 7*d*, 9*b*, 9*d*}, *i* = 1, · · · , *k*, and, *k* = 15.

#### **4. Conclusions**

In this article, various aspects of some parameters of the GA for the attack on block ciphers were studied. In the first place, a way of estimating the time that the GA takes in a given number of generations was proposed, having an average of the time that this algorithm takes in one generation. This study is important to jointly evaluate different parameters and make the best decisions according to the computational capacity, available time, and an adequate selection of the size of the search space when using the BBM and TBB methodologies. On the other hand, several fitness functions were proposed with favorable results in the experiments with respect to the fitness functions using only the Hamming distance. In this sense, it was found that the fitness functions that use the decimal distance, in general, are more efficient than those that use only the Hamming distance, especially in the methodology BBM.

As future work, several directions are possible. Similar studies can be carried out with the GA working with other parameters, such as varying the crossover probability and mutation rate and making comparisons regarding the percentage of success of the method. It is also recommended to explore other heuristic techniques and to evaluate the use of whole space partitioning methods so that the methods work closed on the subsets. In the same way, it is also recommended to investigate the combined use with some other tools such as machine learning, deep learning, ANN, SVM, and GEP.

**Author Contributions:** Conceptualization, O.T.-C. and M.B.-Q.; methodology, O.T.-C., M.B.-Q., M.A.B.-T., and G.S.-G.; software, O.T.-C.; validation, O.T.-C., M.B.-Q., M.A.B.-T., O.R., and G.S.-G.; formal analysis, M.A.B.-T., M.B.-Q., O.R., and G.S.-G.; investigation, O.T.-C., M.B.-Q., M.A.B.-T., and G.S.-G.; writing—original draft preparation, O.T.-C. and M.B.-Q.; writing—review and editing, O.T.-C., M.B.-Q., M.A.B.-T, O.R., and G.S.-G.; visualization, O.T.-C.; supervision, M.A.B.-T., M.B.-Q., O.R., and G.S.-G. All authors read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Informed Consent Statement:** Informed consent was obtained from all subjects involved in the study.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


## *Article* **Sensitivity Analysis of Key Formulations of Topology Optimization on an Example of Cantilever Bending Beam**

**Martin Sotola 1,2 , Pavel Marsalek 1,2,\* , David Rybansky 1,2 , Martin Fusek <sup>1</sup> and Dusan Gabriel <sup>2</sup>**


**Abstract:** Topology optimization is a modern method for optimizing the material distribution in a given space, automatically searching for the ideal design of the product. The method aims to maximize the design performance of the system regarding given conditions. In engineering practice, a given space is first described using the finite element method and, subsequently, density-based method with solid isotropic material with penalty. Then, the final shape is found using a gradientbased method, such as the optimality criteria algorithm. However, obtaining the ideal shape is highly dependent on the correct setting of numerical parameters. This paper focuses on the sensitivity analysis of key formulations of topology optimization using the implementation of mathematical programming techniques in MATLAB software. For the purposes of the study, sensitivity analysis of a simple spatial task—cantilever bending—is performed. This paper aims to present the formulations of the optimization problem—in this case, minimization of compliance. It should be noted that this paper does not present any new mathematical formulas but rather provides an introduction into the mathematical theory (including filtering methods and calculating large-size problems using the symmetry of matrices) as well as a step-by step guideline for the minimization of compliance within the density-based topology optimization and search for an optimal shape. The results can be used for complex commercial applications produced by traditional manufacturing processes or by additive manufacturing methods.

**Keywords:** topology optimization; optimization; filtering; method; penalization; weight factor; FEM; MATLAB; SIMP

#### **1. Introduction**

Topology optimization is a calculation of the distribution of materials within a structure without a known pre-defined shape. This distribution calculation yields a "black and white pattern" where black places indicate full material while white places represent voids (i.e., places where material can be removed). Because the distribution is solved over a general region, topology optimization allows us to acquire a unique, innovative, and effective structure. The principle of topology optimization is presented on the example of cantilever bending in Figure 1, where the initial geometry (given space) is depicted on the left and the optimal shape on the right. As apparent from Figure 1, topology optimization plays nowadays an important role in engineering practice. Usually, it allows the designer to reduce the weight of the part without losing too much of its previous properties such as stiffness, natural frequency, etc.

**Citation:** Sotola, M.; Marsalek, P.; Rybansky, R.; Fusek, M.; Gabriel, D. Sensitivity Analysis of Key Formulations of Topology Optimization on an Example of Cantilever Bending Beam. *Symmetry* **2021**, *13*, 712. https://doi.org/ 10.3390/sym13040712

Academic Editors: Sun Young Cho and Aviv Gibali

Received: 11 March 2021 Accepted: 16 April 2021 Published: 18 April 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

**Figure 1.** Topology optimization of cantilever bending; the initial geometry (**left**) and the optimal shape in the form of design variables (**right**).

Issues associated with topology optimization are studied by many engineers and researchers. The Finite Elements Method (FEM) is the most widely used technique for the analysis of discretized continuum. Femlab [1], FreeFem++ [2,3] and ToPy [4] are rapidly growing engineering tools supporting topology optimization. The problem of topology optimization was described, e.g., by Bendsøe and Sigmund [5–8]. Many computer tools have been prepared, including tools in the MATLAB platform (MathWorks, Natick, MA, United States of America). Liu et al. [9] described a three-dimensional (3D) topology optimization using MATLAB scripts. They described the necessary steps of optimization and provided scripts for individual steps. Their scripts are freely available and can be modified in accordance with the authors' instructions. It should be noted that although the scripts have great educational value, their practical usage is limited as they are applicable only for simple shapes and cannot work with imported meshes. The same can be said about the paper by Sigmund et al. [7] who investigated two-dimensional optimization and introduced sensitivity filtering. Master thesis by William Hunter [10] worked on 3D topology optimization. The author described in depth the theoretical background as well as its implementation in Python.

Even though topology optimization might seem novel, the first mention of structural topology optimization dates back to 1904 [11]. However, major progress has come only in the last 30 years due to the development of computers and the advancement of new technological processes (in particular, additive manufacturing). Hunar et al. [12] and Pagac et al. [13] illustrated the significance of topology optimization for designing a 3D printed part. Currently, however, topology optimization is perceived by most users as a "black box" producing always correct, i.e., optimal, shapes.

This is, however, largely not true and deeper understanding of the parameters and settings is needed to yield optimal results. A thorough introduction to the problem and a complex guideline for performing sensitivity analysis that would help researchers and engineers with determining correct settings is, however, not available in the current literature. For this reason, we decided to provide such a guideline in this paper.

Hence, the presented paper studies the effect of key formulations of the topology optimization problem on the design performance. In addition, it recommends the values of individual numerical parameters using a cantilever beam problem. Even though it is demonstrated on a simple example, this insight can be used for complex problems of engineering practice. For example, our group [14] described topology optimization of a transtibial bed stump using a custom MATLAB script. Performing sensitivity analysis of the key formulations is recommended for supporting the robustness of the computations in any new problem. Failure to perform so such an analysis may lead to the production of "false-optimal" results. The individual steps of this study were performed in MATLAB as well as in ANSYS Workbench 2019 R3 (ANSYS, Inc., Canonsburg, United States of America, AWB), which helped with the evaluation of the results (assessment of the similarity in resulting values or shapes). The preliminary results of the work are presented in the Master's thesis by Sotola [15].

#### **2. Materials and Methods**

The procedure for calculating the optimal shape of the structure can be divided into three stages: (i.) Preparative stage (ii.) Optimization, and (iii.) Postprocessing, see Figure 2. MATLAB scripts were prepared for each of these stages to automatize the procedure.

The first step of the initiation of optimization lies in preparing the finite element analysis. In this step, the boundary conditions and local stiffness matrices of elements are set up. This paper does not describe this stage in depth because it has been already described in many books; for example, Hughes [16] describes the preparation of local and global stiffness matrices. Still, some advanced recommendations are presented in this paper. At this stage, the global stiffness matrix is also assembled and the reference values of the objective function are calculated. Before optimization, it is also necessary to prepare an initial approximation of volume, i.e., the structural elements are assigned a new material model containing individual design variables for each element. Elements are assigned new values of elasticity modulus, which affect the global stiffness matrix and leads to a new value of the objective function with each iteration.

Optimization itself follows, during which new values of design variables are determined. Subsequently, the terminating criterion is queried and if the process is not terminated based on the criterion being met or the maximal number of iterations exceeded, the process is repeated with new design variables.

In the last stage, the results are recalculated and prepared in the *vtk* format, which can be viewed in the open-source software ParaView (Kitware, Inc., Clifton Park, NY, USA). This paper focuses on the first and second stages; we provide the results of the entire process but do not describe the postprocessing in detail.

**Figure 2.** Diagram of the procedure for topology optimization.

#### *2.1. Description of the Optimization Problem*

As mentioned already, being able to describe the problem as a mathematical input is the key to the full understanding of topology optimization. Such understanding is needed for proper definition of the design variables, the objective function, and the constraint function (mostly inequality constraints). The most widely used method for solving multivariate optimization is called "Karush–Kuhn–Tucker conditions" [17] (also known as Kuhn–Tucker conditions or just KKT conditions).

#### 2.1.1. The Optimization Problem

The objective of the optimization presented in this paper is to minimize compliance arising during volume reduction defined as volume fraction vfrac. The volume fraction is calculated as the ratio of the proposed volume and the original volume. The volume fraction ranges from 0% to 100%.

Several methods can be used for solving the problem of minimization of compliance. For example, the homogenization method [18,19] uses microperforated composites as a base material for shape optimization. Since the number of holes within the domain is not

limited, it can be seen as a method of topology optimization. In another approach, the phase-field method, the domain consists of two "phases", the "void" and the fictional "liquid" which interacts with loads [20,21]. One of the newer methods (meaning newly implemented in commercial software) is the Level-set method. Optimization in this method is solved "above" the fixed domain with a fictional velocity [22].

The density-based method is the last of the methods commonly used for solving the topology problem. This method can be viewed as mature and is easy to implement. In this paper, we focused on this method, which uses a continuous design variable ranging from 0% to 100%, see Figure 1. In the literature, the design variable is usually referred to as the density; note, however, that this "density" of the element has no clear physical meaning. In 2D space, the density can be pictured as a variable thickness of sheet metal but in 3D space, it is not easy to assign a tangible meaning to this term. The Solid Isotropic Material with Penalty (SIMP) model is a popular interpolation scheme for definition of material that would be subsequently used in the density-based method. In this model, the elasticity modulus is in a power law relation with the design variable and can be described using the following equation.

$$E\_{\mathfrak{e}} = \widetilde{\mathfrak{x}}\_{\mathfrak{e}}^{p} \cdot E\_{0\prime} \ \widetilde{\mathfrak{x}}\_{\mathfrak{e}} \in < 0, 1 >, \tag{1}$$

where *<sup>E</sup>*<sup>0</sup> is the elastic modulus of the base material, *<sup>p</sup>* is the penalization and *<sup>x</sup>*e*<sup>e</sup>* is the "modified" filtered design variable. The reasons for the "modification" are described in the following Sections 2.1.2 and 2.1.6.

After preparation of the material model, the objective function of minimizing compliance is defined as

$$\min: \mathcal{c}(\widetilde{\mathfrak{x}}\_{\mathfrak{e}}) = \{f\}^T \cdot \{\mathfrak{u}(\widetilde{\mathfrak{x}}\_{\mathfrak{e}})\},\tag{2}$$

where *c* is the deformation energy, { *f* } is the force vector and {*u*} is the displacement. One could argue that the deformation energy of the linear material should be divided by two; however, as it is only a scalar variable, it will not affect the optimization itself. The constraint Equation (or, rather, inequality in this case) is defined as

$$\mathbf{v}(\tilde{\mathbf{x}}\_{\varepsilon}) = \{\,\,\widetilde{\mathbf{x}}\_{\varepsilon}\}^T \cdot \{v\_{\varepsilon}\} - \mathbf{v}\_{\text{frac}} \cdot \sum \{v\_{\varepsilon}\} \le \mathbf{0},\tag{3}$$

where {*ve*} is the vector containing the volume of each element. Displacement is solved from the equilibrium equation

$$\{f\} = [\mathbf{K}] \cdot \{\boldsymbol{u}\},\tag{4}$$

where [*K*] is the global stiffness matrix. Assembly of stiffness matrices is described in Section 2.1.4. The design variable is defined as

$$0 \le \tilde{\mathfrak{X}}\_{\ell} \le 1. \tag{5}$$

Equation (2) can be written in a simpler form after applying SIMP

$$\min: \mathcal{L}(\widetilde{\boldsymbol{x}}\_{\varepsilon}) = \sum\_{\varepsilon=1}^{n} \mathbb{E}\_{0} \, \widetilde{\boldsymbol{x}}\_{\varepsilon}^{p} \left\{ \boldsymbol{u}\_{\varepsilon} \right\}^{T} [\mathbf{k}\_{0}] \{\boldsymbol{u}\_{\varepsilon}\} \tag{6}$$

where {*ue*} is the vector of element displacement, [k0] is the stiffness matrix with unit elastic modulus and *n* is the number of elements. Penalization is introduced in Section 2.1.3.

One can notice that the problem is defined as minimizing the deformation energy, which leads to the higher stiffness of the structure. Using the KKT conditions, the Lagrange multipliers convert the constrained problem to the equivalent unconstrained problem [17,23]. The first step is assembling the Lagrange function

$$L(\widetilde{\mathfrak{x}}\_{\mathfrak{e}}, \lambda) = \mathfrak{c}(\widetilde{\mathfrak{x}}\_{\mathfrak{e}}) + \lambda \cdot \mathfrak{v}(\widetilde{\mathfrak{x}}\_{\mathfrak{e}}),\tag{7}$$

where *<sup>c</sup>*(*x*e*e*) is the objective function, *<sup>v</sup>*(*x*e*e*) is the constraint function and *<sup>λ</sup>* the Lagrange multiplier. It is necessary to solve the derivation of the Lagrange function with respect to the design variables because it is necessary to find a stationary point; at that point, the derivation is equal to 0

$$\nabla L(\widetilde{\mathfrak{x}}\_{\varepsilon}, \lambda) = \frac{\partial c(\widetilde{\mathfrak{x}}\_{\varepsilon})}{\partial \widetilde{\mathfrak{x}}\_{\varepsilon}} + \lambda \cdot \frac{\partial v(\widetilde{\mathfrak{x}}\_{\varepsilon})}{\partial \widetilde{\mathfrak{x}}\_{\varepsilon}} = 0. \tag{8}$$

The following equation is known as the complementary slackness condition, determining whether the constraint function is active or passive

$$
\lambda \cdot v(\tilde{\mathfrak{x}}\_{\varepsilon}) = 0. \tag{9}
$$

The next condition is that the Lagrange multiplier is not negative

$$
\lambda \ge 0.\tag{10}
$$

The last condition is the derivation of the Lagrange function with respect to Lagrange multipliers. In simple terms, it determines whether or not the particular point is the KKT point, i.e., the optimum (there are fundamental theorems proving that the solution is automatically optimal, see more in [17]). This condition is not important for solving the problem, it is important for results evaluation. Rearranging the equation and adding the variable *B<sup>e</sup>* leads to the equation

$$B\_{\mathcal{E}} = 1 = \frac{-\frac{\partial c(\widetilde{\mathcal{X}}\_{\mathcal{E}})}{\partial \widetilde{\mathcal{X}}\_{\mathcal{E}}}}{\lambda \frac{\partial v(\widetilde{\mathcal{X}}\_{\mathcal{E}})}{\partial \widetilde{\mathcal{X}}\_{\mathcal{E}}}}.\tag{11}$$

It is obvious that the optimum of the element is met if *B<sup>e</sup>* = 1. Preparing the first derivation of the objective function determined by Equation (2) with respect to the design variables can be tricky as at the first sight, there is no evident "influence" of the design variables. It can be solved numerically but this would require a calculation of the displacement for every individual possible design variable. However, it is possible to calculate derivation using the adjoint method [5], i.e., to add the equilibrium equation into the objective function

$$\mathcal{L}(\tilde{\mathbf{x}}\_{\varepsilon}) = \{f\}^T \{\boldsymbol{\mu}(\tilde{\mathbf{x}}\_{\varepsilon})\} + \{\boldsymbol{\eta}\}^T ([\![\mathbf{K}]\{\boldsymbol{\mu}(\tilde{\mathbf{x}}\_{\varepsilon})\}\!] - \![\![f]\!])\_{\prime} \tag{12}$$

where *η* is a vector of non-zero variables (also unknown for now). Adding the equilibrium equation does not change the objective function because it is equal to zero. Similarly, the vector *η* does not change the function. Let us assume that the exterior forces are not independent on the design variables. Then, the first-order derivation of the objective function is

$$\frac{\partial c}{\partial \{\widetilde{\boldsymbol{x}}\_{\varepsilon}\}} = \{f\}^{T} \frac{\partial \{\boldsymbol{u}\}}{\partial \{\widetilde{\boldsymbol{x}}\_{\varepsilon}\}} + \{\eta\}^{T} \left(\frac{\partial [\boldsymbol{K}]}{\partial \{\widetilde{\boldsymbol{x}}\_{\varepsilon}\}} \{\boldsymbol{u}\} + [\boldsymbol{K}] \frac{\partial \{\boldsymbol{u}\}}{\partial \{\widetilde{\boldsymbol{x}}\_{\varepsilon}\}}\right). \tag{13}$$

Rearranging the previous Equation (13) leads to factoring out the derivation of the displacement

$$\frac{\partial \mathcal{C}}{\partial \{\widetilde{\boldsymbol{x}}\_{\varepsilon}\}} = \left( \{f\}^{T} + \{\boldsymbol{\eta}\}^{T}[\boldsymbol{K}] \right) \frac{\partial \{\boldsymbol{u}\}}{\partial \{\widetilde{\boldsymbol{x}}\_{\varepsilon}\}} + \{\boldsymbol{\eta}\}^{T} \frac{\partial[\boldsymbol{K}]}{\partial \{\widetilde{\boldsymbol{x}}\_{\varepsilon}\}} \{\boldsymbol{u}\}. \tag{14}$$

To get rid of the derivation of displacement, it is critical that the term in brackets (adjoint equation) is equal to zero. This means that our added unknown variable has to be equal to

$$\{\eta\}^T = -\{f\}^T[\mathbb{K}]^{-1} = -\{\boldsymbol{\mu}\}^T.\tag{15}$$

It is apparent that the vector {*η*} is already solved. The final form of the derivation of the objective function (with added SIMP model) is

$$\frac{\partial c(\widetilde{\mathbf{x}}\_{\varepsilon})}{\partial \widetilde{\mathbf{x}}\_{\varepsilon}} = -E\_0 \, p \, \widetilde{\mathbf{x}}\_{\varepsilon}^{(p-1)} \, \{\boldsymbol{u}\_{\varepsilon}\}^T [\mathbf{k}\_0] \{\boldsymbol{u}\_{\varepsilon}\}. \tag{16}$$

The first order derivation of the constraint function is defined asIt should be noted that in the literature (for example, [8]), it is possible to find another equation describing the constraint function and its derivation

$$\frac{\partial v(\tilde{\mathbf{x}}\_{\varepsilon})}{\partial \tilde{\mathbf{x}}\_{\varepsilon}} = \{v\_{\varepsilon}\}. \tag{17}$$

It should be noted that in the literature (for example, [8]), it is possible to find another equation describing the constraint function and its derivation

$$v(\tilde{\mathbf{x}}\_{\varepsilon}) = \frac{\{\tilde{\mathbf{x}}\_{\varepsilon}\}^T \cdot \{v\_{\varepsilon}\}}{\sum \{v\_{\varepsilon}\}} - \mathbf{v\_{\text{frac}}} \le \mathbf{0},\tag{18}$$

$$\frac{\partial v(\tilde{\mathbf{x}}\_{\varepsilon})}{\partial \tilde{\mathbf{x}}\_{\varepsilon}} = \frac{\{v\_{\varepsilon}\}}{\sum \{v\_{\varepsilon}\}}.\tag{19}$$

That form of the equation usually depends on the solver and its settings. For example, "MMA-based" solvers prefer the constraint function defined by Equations (18).

#### 2.1.2. Material Model

In this paper, the density-based method used the Solid Isotropic Material with Penalty (SIMP) material model, which is a power-law relation of the design variables. The SIMP material model is used in solvers such as ANSYS (ANSYS, Inc., Canonsburg, PA, USA), MSC Nastran (MSC Software, Irvine, CA, USA), etc. The elastic modulus of the element is defined as

$$E\_{\mathfrak{e}} = \mathfrak{x}\_{\mathfrak{e}}^p \cdot E\_{0\prime} \ \mathfrak{x}\_{\mathfrak{e}} \in < \mathbf{0}, 1 >, \tag{20}$$

where *E*<sup>0</sup> is the elastic modulus of the base material, *p* is the penalization and *x<sup>e</sup>* is the "unmodified" design variable. This definition looks very similar to Equation (1); however, here, we use unfiltered designed variable to give a clearer picture of the need for filtering (see below).

During topology optimization, many problems can occur. For example, the so-called checkerboard pattern problem [24–26] is very common. This title describes the distribution of the structural elements in a checkerboard-like arrangement in certain areas of the part. Figure 3 shows the checkerboard patterns on the cantilever beam.

**Figure 3.** Checkerboard pattern problem (**left**) and its possible solution (**right**).

In our case, i.e., when discussing the problem of minimizing compliance, the computation may consider such a solution to be ideal as the checkerboard pattern creates artificial regions with higher stiffness. However, as obvious from Figure 3, this is not true and, moreover, the part cannot be manufactured. Another common problem is the insufficient number of connected nodes (four or more nodes are needed for the hexahedral element), which causes the formation of possible joints in the structure and, again, makes it impossible to manufacture.

The mesh-dependency of optimization results is a crucial problem [5,26]. As the name suggests, this problem results from the used discretization and its refinement. In context of

structure stiffness, the reason is quite simple—increasing the number holes in the structure without changing its volume leads to increase of stiffness. Finer meshes facilitate this operation—they allow us to create a higher numbers of holes and, therefore, are capable of providing different (superior) results than coarser meshes. On the other hand, finer meshes might result in more complex structures that are difficult to manufacture.

The last important problem of topology optimization is the non-uniqueness of solutions, which makes results evaluation trickier [27]. Generally, the problem of optimization can have a single solution or up to an infinite number of solutions depending on whether the problem is convex or not.

One of the ways of solving the above-mentioned problems is to use a suitable filter(s). This solution has not been mathematically confirmed, but many numerical experiments suggest that the results could be considered optimal [7]. The description of each of the filters used in our study is presented in Section 2.1.6.

To prevent numerical difficulties, the modulus of void (passive) elements is introduced into the material model. This modification helps to reduce the risk of having a singular stiffness matrix. The final equation of the material model is

$$\begin{aligned} E\_{\varepsilon} &= E\_{\min} + \widetilde{\mathbf{x}}\_{\varepsilon}^{p} \cdot (E\_0 - E\_{\min}) \\ \widetilde{\mathbf{x}}\_{\varepsilon} &= \widetilde{\mathbf{x}}\_{\varepsilon}(\mathbf{x}\_{\varepsilon}) \,, \ \mathbf{x}\_{\varepsilon} \in < \mathbf{0} \, \mathbf{1} > \,, \end{aligned} \tag{21}$$

where *<sup>E</sup>min* is the elasticity modulus of passive elements and *<sup>x</sup>*e*<sup>e</sup>* is the filtered design variable (density).

#### 2.1.3. Penalization

The numerical scheme should lead to a black & white design (or 1-0 design with white to be removed). One of the possible approaches is to ignore the physical importance of elements with intermediate density (grey areas) and consider them "black", leading to their preservation. However, the physical relevance is discussed a lot since many interpolation methods can remove further parts of the grey regions. If the optimization is prematurely terminated, the stiffness (compliance) of the grey areas plays an important role in the evaluation of results. This issue is discussed by Bendsøe [5,6]. As mentioned above, the SIMP material model is suitable for FEM optimization as it assigns elasticity to each element. However, it can be used as the material model only if the penalization meets the following criteria

$$p \ge \max\left\{\frac{2}{1-\mu}, \frac{4}{1+\mu}\right\}, \text{ (in 2D)}; \ p \ge \max\left\{15\frac{1-\mu}{7-5\mu}, \frac{3}{2}\frac{1-\mu}{1-2\mu}\right\}, \text{ (in 3D)}\tag{22}$$

where *µ* is the Poisson's ratio. This means that different values of penalization must be calculated for each Poisson's ratio. The resulting values of the penalization for volume elements with the following Poisson's ratio are

$$p\left(\mu = \frac{1}{5}\right) \ge 2, \ p\left(\mu = \frac{1}{3}\right) \ge 3, \ p\left(\mu = \frac{2}{5}\right) \ge 4.5. \tag{23}$$

In this paper, the base material has a Poisson ratio of *µ* = 0.30; therefore, the penalization *p* = 3 was chosen for the following topology optimization.

#### 2.1.4. Finite Element Method

From our experience, it is recommended for MATLAB implementation that the stiffness matrices is prepared with unit elastic modulus and saved in memory. The stiffness matrix of the solid element is volume integrated using the stress–strain matrix of the material [*C*0] with the unit elasticity modulus

$$\mathbb{E}\left[\mathbf{k}\_0\right] = \iiint \left[\mathbf{B}\right]^\mathrm{T} \left[\mathbf{C}\_0\right] \left[\mathbf{B}\right] \,dV\_\prime \,d\mathbf{V}\_\prime \tag{24}$$

where [*B*] is the strain–displacement matrix [28]. The integral is usually solved numerically. Before the assembly of the global stiffness matrix, every local stiffness matrix is multiplied by the corresponding elasticity modulus (i.e., the modulus of the SIMP material model).

$$[\mathbf{k}\_{\mathfrak{e}}(\tilde{\mathbf{x}}\_{\mathfrak{e}})] = E\_{\mathfrak{e}}(\tilde{\mathbf{x}}\_{\mathfrak{e}}) \cdot [\mathbf{k}\_{\mathbf{0}}].\tag{25}$$

The MATLAB implementation of calculating the stiffness matrix of linear elements is presented by Bhatti [29,30]. With some effort, the procedures were modified to the calculation of quadratic elements. The procedures use full integration of the individual elements using Gauss points (2 × 2 × 2 scheme for linear elements, 3 × 3 × 3 for quadratic elements).

For testing purposes, meshes were prepared using the AWB software [31]. They were constructed either of tetrahedral (TET) or hexahedral (HEX) elements (i.e., no result with mixed mesh was evaluated in this paper). The calculation of the stiffness matrices of linear elements is fast thanks to the low number of degrees of freedom (DOF); in effect, solving a single load case was swift when using this approach; however, it comes with a disadvantage in the form of locking as linear elements with full integration have a tendency to shear and volume locking [28]. There are numerous ways of fixing this problem; nonetheless, in this paper, the shear locking effects are "neglected" during the optimization but are mentioned in results. Meshes are separated into three groups: coarse, normal, and fine meshes according to the size (length) of individual elements.

#### 2.1.5. The Optimality Criterion Algorithm

The design variable of the elements was updated using an algorithm called the Optimality criterion (OC) [5,32]. The name indirectly refers to the used method, i.e., KKT conditions. The algorithm is also implemented in the AWB software. To find the material distribution of the structure, a fixed updating scheme is proposed as

$$\mathbf{x}\_{\varepsilon}^{\text{new}} = \begin{cases} \max(0, \mathbf{x}\_{\varepsilon} - m), & \text{if } \mathbf{x}\_{\varepsilon} \cdot \mathbf{B}\_{\varepsilon}^{\eta} \le \max(0, \mathbf{x}\_{\varepsilon} - m) \\ \min(1, \mathbf{x}\_{\varepsilon} + m), & \text{if } \mathbf{x}\_{\varepsilon} \cdot \mathbf{B}\_{\varepsilon}^{\eta} \ge \min(1, \mathbf{x}\_{\varepsilon} + m) \\ \mathbf{x}\_{\varepsilon} \cdot \mathbf{B}\_{\varepsilon}^{\eta} & \text{otherwise} \end{cases} \tag{26}$$

where *B<sup>e</sup>* is constructed using Equation (11). As already mentioned above, the optimum of the element is found if *B<sup>e</sup>* = 1. In other words, the design variable increases if *B<sup>e</sup>* > 1 and decreases if *B<sup>e</sup>* < 1. Changing the move limit *m* and tuning parameter *η* can lead to a lower number of iterations. Bendsøe [5] recommended values of *η* = 0.5 and *m* = 0.2. New values of design variables depend on the Lagrange multiplier, which has to be solved in the inner loop to ensure that the constraint function is satisfied. This leads to a reduction of the multivariate problem to one-dimensional (1D) optimization, which can be solved by various methods, such as the bisection method, golden section search, or methods using derivation (such as the Newton–Raphson method or secant method) [17,23]. In this work, the bisection method obtained from the paper by Liu [9] was used as the 1D optimization method.

A detailed description of the Optimality criterion method is presented in the dissertation thesis by Munro [32] and the Master's thesis by Hunter [4]. These theses describe the relationships between the Optimality criterion and Sequential approximate optimization (SAO). The SAO method can be solved using the duality principle. The purpose of this method is to find an equivalent subproblem (dual problem) that is easier to solve than the primary problem. The method is used for preparing a scheme (similar to the OC fixed scheme) with only one constraint function, which can be subsequently solved. The primary problem (in this case, minimizing compliance) is then reduced to maximizing the Lagrange multiplier.

2.1.6. Filtering Methods

Density filtering is one of the methods of solving the above-mentioned problems (such as the checkerboard pattern). A common parameter of the filters, weight factor *Hij*, is defined as

$$H\_{ij} = \begin{cases} \mathbb{R} - dist(i, j) & \text{if } dist(i, j) \le \mathbb{R} \\ 0 & \text{if } dist(i, j) > \mathbb{R} \end{cases} \tag{27}$$

where *R* is the radius of the filter, *dist*(*i*, *j*) is the operator calculating the distance between the center of an element *i* and the center of an element *j*. This type of weight factor is called linear. If the distance between elements is greater than the radius, the weight factor is equal to zero; if the distance is equal to zero, the weight factor is equal to the size of the radius. In Figure 4, 2D examples are presented; in 2D, the radius defines a circular neighborhood. In a 3D problem, the radius forms a sphere.

**Figure 4.** 2D Neighborhood of element (**left**), examples of the radius dependence on the element size *ES* in 2D; R= 1.2 *ES* for the red circle, R= 1.5 *ES* for the blue circle, R= 2.0 *ES* for the yellow circle and R= 3.0 *ES* for the purple circle.

An alternative approach is to use the normal distribution (Gauss function) [33,34]. Compared to the linear function, the Gauss function is smoother but in reality, there might not be a real benefit in using this alternative [8]. The weight factor is defined as

$$H\_{ij} = e^{-\frac{1}{2} \left( \frac{3 \cdot \text{dist}(i,j)}{R} \right)^2}. \tag{28}$$

In Figure 5, both functions are displayed. Due to the different characters of values of the weight factor, the linear function had to be normalized by the radius to have values ranging from 0 to 1, the same as the Gauss function. Furthermore, the distance was normalized to the radius (i.e., divided by the radius to ensure independence on it).

Preparing a matrix of weight factors can be challenging if the mesh is imported (if the mesh is not imported but created by additional scripts, the weight factor can be calculated during the mesh preparation). The following should be pointed out: Firstly, the preparation should use a function for creating sparse matrices because the matrix of weight factors is mainly sparse. Secondly, if the distance between the *i*-th and *j*-th element is constant, there is no need to create a "full" matrix but only the upper triangle of the matrix [*Htri*]. Non-diagonal components of the matrix [*Htri*] are multiplied by 2. The following equation uses "symmetry" for composing the full matrix

$$\left[H\right] = \frac{1}{2} \cdot \left( \left[H\_{tri}\right] + \left[H\_{tri}\right]^T \right) \tag{29}$$

where *T* is the operator of the transpose. Lastly, with the growing amount of elements, the process of preparation is becoming more time-consuming even though calculating

only "half" of the matrix. This basic approach is, up to 10 k elements, fast. However, the time of solution is growing exponentially for meshes with over 10k elements. In such cases, it is recommended to invest time into finding an appropriate method for speeding up the preparation. This could be done by dividing the mesh into mutually overlapping subzones with less than 10 k elements to ensure fast calculation. These subzones should be solved individually using the parallel toolbox (package). Using this approach, it should be possible to prepare the weight factors of complex meshes in a reasonable time.

Three filters are analyzed in this paper: (i) Density filter, (ii) Sensitivity filter, and (iii) Greyscale filter (G). The Density filter (D) and Sensitivity filter (S) use weight factors while the greyscale filter is an addition to the OC scheme.

#### Density Filter

The fundamental function of density filtering is

$$\widetilde{\mathbf{x}}\_{i} = \frac{\sum\_{j=1}^{n} H\_{ij} \mathbf{x}\_{j}}{\sum\_{j=1}^{n} H\_{ij}}, \; i = 1, \dots, n,\tag{30}$$

where *<sup>x</sup>*e*<sup>i</sup>* is the filtered design variable (density) of the *i*-th element, *Hij* is the weight factor and *n* is the number of elements [35–37]. Every equation containing the design variables has to be adjusted to allow filtering (including partial derivations)

$$\frac{\widehat{\partial c}}{\partial \mathbf{x}\_{\varepsilon}} = \frac{\sum\_{j=1}^{n} H\_{i\bar{j}} \frac{\partial c}{\partial \mathbf{x}\_{\varepsilon}}}{\sum\_{j=1}^{n} H\_{i\bar{j}}}, \frac{\widehat{\partial v}}{\partial \mathbf{x}\_{\varepsilon}} = \frac{\sum\_{j=1}^{n} H\_{i\bar{j}} \frac{\partial v}{\partial \mathbf{x}\_{\varepsilon}}}{\sum\_{j=1}^{n} H\_{i\bar{j}}}. \tag{31}$$

In this case, the design variable of the element is averaged over its neighborhood. It ensures the smoothing of stiffness. The filtered density is applied during the construction of the global stiffness matrix before solving the static structural analysis. Generally, the density tends to have a value of 0 or 1 (which could be ideal), but after applying the filter, regions of intermediate (grey density) design variable appear, which are then penalized by the SIMP model. Deformation energy is also averaged and shared with the neighborhood of elements.

#### Sensitivity Filter

This filter is based on filtering of the sensitivity (i.e., of the first partial derivation of Lagrange function). Experience has proven that filtering sensitivity ensures mesh independence of results and is time-effective [7]. It is also easy to implement and does not increase the complexity of the problem. The filter is purely heuristic but has been proven to yield similar results as the gradient constraint method [5]. The fundamental equation of this filter is

$$\frac{\widehat{\partial c}}{\partial \mathbf{x}\_{i}} = \frac{1}{\max(10^{-3}, \mathbf{x}\_{i}) \cdot \sum\_{j=1}^{n} H\_{ij}} \cdot \sum\_{j=1}^{n} H\_{ij} \cdot \mathbf{x}\_{i} \cdot \frac{\partial c(\mathbf{x}\_{i})}{\partial \mathbf{x}\_{c}}.\tag{32}$$

Grayscale Filter

The last filter is an addition to the previous filters called the Grayscale filter [38]. Applying this filter should reduce the grey areas. The objective of the filter is also to penalize the volume constraint in the OC algorithm. The parameter *q* is implemented in the OC scheme and its value can be constant or gradually increasing by multiplication (in our case, coefficient *q* was multiplied by 1.01 each iteration). The amended OC scheme is described as

$$\mathbf{x}\_{\varepsilon}^{\mathbf{new}} = \begin{cases} \max(0, \mathbf{x}\_{\varepsilon} - m), & \text{if } \mathbf{x}\_{\varepsilon} \cdot \mathbf{B}\_{\varepsilon}^{\eta} \le \max(0, \mathbf{x}\_{\varepsilon} - m) \\ \min(1, \mathbf{x}\_{\varepsilon} + m), & \text{if } \mathbf{x}\_{\varepsilon} \cdot \mathbf{B}\_{\varepsilon}^{\eta} \ge \min(1, \mathbf{x}\_{\varepsilon} + m) \\ (\mathbf{x}\_{\varepsilon} \cdot \mathbf{B}\_{\varepsilon}^{\eta})^{q} & \text{otherwise} \end{cases} \tag{33}$$

This filter is activated after 15 iterations of optimization. Usually, it is limited by the maximum value of the coefficient *q* (in this paper, the maximum was set to 5). If the q coefficient maximum was set to 1, the filter would be deactivated. The idea of the filter is to underestimate the intermediate density leading to a zero value (void density). Underestimation occurs also in the inner iteration. This leads to fewer iterations needed for finding solutions.

#### 2.1.7. Displacement Solver and Termination Criteria

To solve the displacement Equation (4), it is possible to use the direct solver present in MATLAB. However, if the number of degrees of freedom is high (above sixty thousand), the pursuit to provide an accurate result is too ambitious. In such cases, therefore, it is appropriate to switch to the iteration solver. MATLAB includes a solver using the conjugate gradient method with preconditioning [39]. It is recommended to use the simplest preconditioning matrix, the diagonal (Jacobi) preconditioner. The reason is that to ensure the best possible solving stability, assembly of the preconditioner is needed in every iteration and, hence, the simpler is the preconditioner, the faster is the solution. Use of a different (more complex) preconditioner, such as Incomplete Cholesky factorization (with various settings),could reduce the number of iterations but the computing time would be the same or higher due to a slow preconditioner assembly (we tested these in the preliminary stage but the detailed results are not presented here). An alternative approach would be to prepare a "universal" preconditioner from the reference matrix only once and to use in every iteration [40]. In the presented study, however, we used the first approach with assembling the diagonal preconditioner each iteration.

It is important to determine the termination criteria. The maximal number of iterations is the first common criterion. Usually, the value ranges between 200 and 500. It should be noted that a complex task (such as a complex geometry or complex loads) requires more iterations but from our experience, 100 iterations are enough for a simple task. The second criterion is defined as the tolerance of sufficient optimization at which the calculation is terminated. There are two possible options.

The first option is to calculate the change in the values of design variables between the current and previous iteration [7]. This change is compared with the chosen tolerance value and once the change in the design variables is below the tolerance value, the calculation is terminated. Thus, if the tolerance is set to a low value, the number of iterations will increase. This is defined as

$$\max \left| \mathbf{x}\_{\mathcal{e}\_i} - \mathbf{x}\_{\mathcal{e}\_{i-1}} \right| \le \varepsilon \tag{34}$$

where *xe<sup>i</sup>* are the design variables of the current iteration, *<sup>x</sup>ei*−<sup>1</sup> are the design variables of the previous iteration and *ε* is the tolerance value.

The other option is to calculate the ratio of the change of the objective function to the current objective function value. It is assumed that the changes of the objective function near the stationary point are minimal. Compared to the first option, the number of iterations is lower (it may be reduced by as much as half). This option is defined as

$$\left| \frac{c(\mathfrak{x}\_{\varepsilon\_i}) - c(\mathfrak{x}\_{\varepsilon\_{i-1}})}{c(\mathfrak{x}\_{\varepsilon\_i})} \right| \le \varepsilon\_\prime \tag{35}$$

where *c*(*xe<sup>i</sup>* ) is the value of the objective function in the *<sup>i</sup>*-th iteration, *<sup>c</sup>*(*xei*−<sup>1</sup> ) is the value of the objective function of the previous iteration (*i*−1,) and *ε* is the tolerance value. This method, however, comes with a risk of premature termination; for this reason, the first option of terminating criterion was used.

#### *2.2. Key Formulations*

It is apparent from the previous sections that the optimization comes with many formulations and parameters. Each parameter can be changed, which might reduce the number of iterations, improve the values of objective functions or of the desired volume fraction; on the other hand, the changes may also lead to solver instability, premature termination, or ineffective shape of the part. In this paper, five key formulations are analyzed and evaluated:

**Formulation of the filter radius** mentioned in Section 2.1.6 is important because the radius defines the element's neighborhood. If the defined range is too small, the energy is distributed only to a few elements. However, if the neighborhood is too large, the energy is scattered to a point where it is difficult to evaluate the optimum.

**Formulation of the filter type** was already mentioned in Section 2.1.6 but the theory does not provide an answer to the question of which filter should perform best.

**Formulation of the penalization** was also mentioned in Section 2.1.3; the theory, however, is able to provide only the lower boundary, not the upper one.

**Formulation of the element approximation** mentioned in Section 2.1.4 is a necessary step in the initiation of optimization. The element approximations greatly affect the accuracy of the solved displacement and the value of the objective function.

**The formulation of the type of the weight factor** mentioned in Section 2.1.6 is defined by two functions. However, theory does not provide enough evidence to decide, which function offers the better performance.

These key formulations were tested on several numerical examples including planar and spatial problems (see Figure 6) but to ensure the clarity of this paper, only one example is presented.

**Figure 6.** Results of numerical examples of topology optimization; from left to right-four point bending, L problem, two-loadcase cantilever plate, beam with square cross-section subjected to torsion.

#### *2.3. Description of Numerical Test*

The sensitivity analysis was performed using a numerical test-cantilever beam. It is a standard problem of mechanics, well-known to every designer. As the optimal shape can be found intuitively, results evaluation is easier.

The boundary conditions were simple. The beam was fixed on one end and the force acted on the other end's edge (bottom edge). The beam had a rectangular cross-section and was made of steel. It was assumed that the forces caused only a small displacement and the original material model was linear isotropic. All finite element meshes were made of solid elements, see Table 1. The authors used linear HEX elements and linear TET elements. In addition, quadratic HEX and TET elements have been calculated and results (mesh statistics) presented in Section 3.4. The geometry and discretization are shown in Figure 7. The figure also contains material and force parameters. The objective of optimization was to minimize the compliance (i.e., to maximize stiffness). In our paper, the constraint was defined as the volume fraction of 30%.

**Figure 7.** Geometry, discretization and parameters of the numerical test-cantilever beam.


#### **3. Results**

Optimization using MATLAB software was performed on multiple meshes. The maximal number of iterations was set to 200. The maximal change of the design variable (i.e., density) was chosen as the terminating criterion. The tolerance was set to 0.01. The cut-off limit of the design variable for element deactivation (white in the Figures) was 0.01 unless stated otherwise. The penalty was set to a constant value of *p* = 3 unless stated otherwise. Unless stated otherwise, the Density filter was used throughout the paper. In tables, two variants of the objective function are shown. The non-normed value indicates the deformation energy. The normed value is the objective function divided by the reference value of the initial objective function (i.e., the objective function describing the original structure before optimization). In other words, the normed value shows how many times the resulting structure is more compliant than the original reference. The normed value in the linear static analysis should be the same (or, at least, similar) regardless of the applied force. Data and shapes presented in this chapter were prepared in Ansys Workbench (AWB).

#### *3.1. Formulation of Filter Radius*

The radius of the filter is an important parameter since it defines the element's neighborhood. In this case, the radius is dependent on the element size (*ES*). For hexahedral

elements, multipliers were set to 1.2 *ES*, 1.5 *ES*, 2.0 *ES*, and 3.0 *ES*, respectively. For clarification, Figure 4 displays the mentioned radiuses. For tetrahedral elements, multipliers were twice as high (i.e., 2.4 *ES*, 3.0 *ES*, 4.0 *ES* and 6.0 *ES*) to prevent possible inactivation of filters.

Results for the hexahedral and tetrahedral meshes are shown in Figure 8, which is displayed in the dominant (planar) view.

**Figure 8.** Final shapes for linear HEX elements (**left**) and linear TET elements (**right**) with various radii, Density filter (D).

From Table 2, it was apparent that the radius heavily affected the distortion of design variables over the design region. This means if the radius was increased, the value of the objective function (both deformation energy and normed value) also increased. Besides, the volume fraction increased due to the distortion. A fine mesh with small radii yielded a great stiffness-volume ratio, the optimized structure was approximately 2.5 times more compliant than the original structure but the volume reduction was as high as 53%.

A few notes: Radius should never be smaller than 1.5 *ES*/2.4 *ES*. A radius such as 1.2 *ES* greatly limited the capability of filters. In the case of uniform mesh, the radius should be appropriately chosen from the range between 1.5 *ES* and 3.0 *ES*. In the case of the non-uniform (tetrahedral) mesh, the radius should be within the range of 2.4 *ES* to 4.0 *ES*. Radii exceeding the upper value of the mentioned ranges make the structure more compliant. In the case of a coarser TET mesh, the radius *R* = 6.0 *ES* caused an over two-fold increase in the deformation energy than radius *R* = 2.4 *ES*. Besides, the shapes of the coarser meshes with larger radii were not acceptable due to the high representation of "grey" areas (see Figure 8).


**Table 2.** Results of optimization with different radius, for meshes with linear HEX elements (first multiplier) or linear TET elements (latter multiplier).

#### *3.2. Formulation of Filter Type*

The importance of filters was already mentioned in the previous section. The following effects were evaluated: the number of iterations, objective function, and volume fraction. This particular aim of our study was to find the appropriate filter leading to a 0–1 design (and low objective function). Five variants were calculated: No filter (NoF), Density filter (D), Sensitivity filter (S), Grayscale filter with Density filter (G), and Grayscale filter with sensitivity filter (SG). Figure 9 shows each filter.

**Figure 9.** Results of individual filtering algorithms: Density filter (D), Sensitivity filter (S), Grayscale filter (G),Combination of the Sensitivity and Grayscale filters (SG).

At the first look, it should be apparent that the SG combination leads in this case to a purely 0-1 design (with the required volume fraction).

Two variants of the radius were selected: In the first variant, the radius was considered to depend on the element size, namely, it was defined as 1.5 *ES*. Results of this approach are displayed in Figure 10 and Table 3. The way of how the filters acted when the radius was not constant over different meshes should be noted.

**Figure 10.** The final shape for linear HEX elements (**left**) and linear TET elements (**right**) with various filters, HEX radius *R* = 1.5 *ES* and TET radius *R* = 3.0 *ES*.

**Table 3.** Results of optimization for various filters for each mesh, dependent radius, filter radius of linear HEX mesh is *R* = 1.5 *ES* an filter radius of linear TET mesh is *R* = 3.0 *ES*.


In the second variant, the radius was not considered to be dependent on the element size and was assigned a constant value of *R* = 4 mm. This radius was defined as a 1.5 multiple of the element size from the coarse tetrahedral mesh. Results displayed in Figure 11 and Table 4 demonstrate the independence of the results on the mesh (results have

a similar shape and value of the objective function). In this variant, TET meshes performed poorly, especially in combination with the Sensitivity filter, which acted unpredictably at best. A coarser TET mesh could result in an acceptable shape (i.e., a shape similar to that derived using the HEX mesh) but finer TET meshes did not lead to volume reduction, but rather to its increase. In addition, the combination of TET mesh with Sensitivity filter often resulted in premature termination of the computation.

**Figure 11.** The final shape for linear HEX elements (**left**) and linear TET elements (**right**) with various filters, radius *R* = 4 mm.


**Table 4.** Results of optimization for various filters for each mesh, independent radius, filter radius *R* = 4 mm.

It should be noted that unwanted effects, such as shear locking, can occur while solving the static structural analysis and cause an increase in the relative error of the objective function. In case of the radius of *R* = 4 mm, the deformation energy of the Density filter was slightly higher for HEX elements than for TET elements. Hence, relative errors were small enough to justify acceptation of the results.

Our results indicate that all filters discussed in this paper are suitable for use with a uniform mesh. The recommended filter combines the Density or Sensitivity filter with the

Greyscale filter to ensure a low number of iterations. In the case of a fine mesh, the Density filter without the Greyscale filter reached a maximal number of iterations while with the Greyscale filter, only 34 iterations were needed. The deformation energies of both results were similar. In the case of a non-uniform mesh, only the Density and Greyscale filters can be used effectively.

One could argue that the combination of the Sensitivity and Greyscale filters got us a perfect black and white design with the required volume fraction (volume reduction of 70%) while increasing the compliance of the structure only approximately 2.1 times; however, the shape was likely to be prone to buckling since compared to other results, the parts were thin.

#### *3.3. Formulation of Penalization*

Penalization is the parameter of the SIMP material model. Its correct choice is crucial as incorrect penalization would invalidate the model. For the Poisson ratio of *µ* = <sup>1</sup> 3 , it is recommended to use the penalization of *p* = 3 and higher as mentioned above.

Figure 12 demonstrates that the penalizations *p* = 1 and *p* = 2 should not be used since the shapes are not fully optimized (due to the premature termination of optimization). Literature, however, does not set the upper limit of this inequality. The values of deformation energy grew more or less predictably up to a penalization value of *p* > 6, see the values in Table 5. With high penalization values, such as *p* > 6, it was, nevertheless, clear that the results were becoming mesh-dependent. This could be caused by the convergence of the solution to a local minimum rather than a global minimum.

One could choose the continuation method with a gradual increase in the penalization value. This approach should help in acquisition of a reasonable solution. It should be, however, mentioned that the continuation method might not be capable of yielding a true "black and white" design, as reported by Stolpe and Svanberg [41]. This paper does not fully study this strategy; the risks are that a too fast increase of the penalization could lead to numerical difficulties that only increase the number of iterations.

**Figure 12.** The final shapes for different values of penalization *p*, Density filter, radius of filter *R* = 2 mm.


**Table 5.** Values of objection function, deformation energy *c* [mJ] for each values of penalization.

#### *3.4. Formulation of Element Approximation*

During the mesh preparation, it is necessary to choose an element approximation (usually linear or quadratic displacement approximation). Literature suggests that problems such as the checkerboard patterns should be less common with quadratic elements [26]. The advantage is that quadratic elements are less stiff than linear ones. However, the need to solve a larger number of unknowns (DOF) is a considerable disadvantage of this approach. Table 6 lists the mesh statistics.



In our case, the optimization settings were slightly altered. Only the Density filter with a radius of *R* = 1.5 *ES* for HEX elements and *R* = 3.0 *ES* for TET. Element size of *ES* = 1.0 mm was used. The maximal number of iterations was set to 100.

Figure 13 shows that the previous statement about the checkerboard pattern being less common with quadratic elements is only partially true. In the case of quadratic HEX elements, a checkerboard pattern was still present in the front part of the structure (although less than when linear elements were used). In the case of TET elements, not even higher order approximations did reduce the checkerboard pattern. Hence, filtering is highly recommended regardless of whether quadratic or linear elements are used. Table 7 obviates that the linear elements are stiffer, i.e., that they provide lower values of the objective function. However, the time needed to solve the problem with quadratic elements is up to thirty times longer than with linear elements.

1.0


**Table 7.** Results of optimization for linear elements and quadratic elements.

#### *3.5. Formulation of Type of Weight Factor*

In the theoretical part, two types of weight factors are mentioned—one characterized by a linear function, the other by the Gauss function. Nevertheless, according to the literature, there is no evidence that the Gauss function offers any advantages over the linear one [8].

In the investigation of this formulation, only two meshes were tested. Both meshes used linear elements with element sizes of *ES* = 1.0 mm and *ES* = 0.5 mm. Two radii were chosen as *R* = 1.5 *ES* and *R* = 2.0 *ES*, respectively. The maximal number of iterations was 100.

The values of the objective function and volume fraction detailed in Table 8, as well as results shown in Figure 14, indicate a significant resemblance between values acquired using Gaussian and linear weight factors. In the case of a normal mesh (*ES* = 1.0 mm), setting the radius of the linear function to *R* = 1.5 *Es* provided similar results as in the case of the Gauss function with a radius of *R* = 2.0 *ES* (see the highlighted values in Table 8). In the fine mesh (*ES* = 0.5 mm), the similarities were more pronounced than in the normal mesh. This means that the Gauss function does not offer any significant advantage over the linear function as their results are very similar and shapes similar to those produced by the Gauss function can be obtained by simply changing the radius in the linear function/solution. As the preparation of mathematical apparatus is simpler with linear function, we decided to prefer linear solution over the one with the Gauss function.


**Table 8.** Results of optimization for different weight factors.

**Figure 14.** The final shape with different weight functions, mesh with element size *ES* = 1.0 mm (**left**), mesh with element size *ES* = 0.5 mm (**right**), Density filter.

#### **4. Discussion**

To fully understand optimization, one should possess advanced experience in math (namely, calculus and finite element method) and computer science (for example, scripting for iteration solvers). One should also be capable of preparing a correct mathematical formulation of the problem, i.e., of determining the objective function of the problem, constraints, and design variables, and of finding equivalent problems, potentially offering an easier solution than the original one.

Five key formulations were analyzed in this paper. The results make it clear that each formulation affects the resulting shape, number of iterations, volume of the part, etc.

The formulation of the filter radius is crucial for determining the element neighborhood. Choosing too small radius might lead to the deactivation of the filter. Choosing a too large radius leads to dissipation of energy in large areas, resulting in too much "grey". Choosing the filter itself can be difficult since the theory does not provide enough knowledge from the application point of view. The SIMP model was used as a simple material model in this study. The penalty value is a key parameter of this model. However, the theory provides only the bottom boundary but does not inform about the upper one. Failure to limit the upper boundary could lead to invalid results that would be highly mesh-dependent. The maximum reasonable penalty value for the steel cantilever in our study was *p* = 6.

The theory also recommends a quadratic approximation of the displacement of the elements. However, it does not provide clear reasons behind this recommendation. Choosing the right approximation might lead to a faster calculation. Lastly, neither the theory nor practice (as demonstrated by our calculations) provide enough evidence or reason for using the Gauss function as the weight function.

It should be noted that there are additional formulations affecting the final shape. For example, various algorithms, described by Zuo et al., can be used in topology optimization [42]. Another approach to this problem could lie in reducing the computing efforts, as reported by Amir et al. [40]. For larger problems, it would be better to prepare a better displacement solver (for example, a parallel displacement solver as suggested, e.g., by Makropoulos et al. [43]. In our paper, however, we did not study other types of structural analysis, such as the heat transfer of flow optimization.

The algorithm used in our scripts, the Optimality Criterion (OC), was designed only for minimizing the compliance and volume constraint. The advantages of the OC algorithm include its simplicity and rapid updating of the design variables. The disadvantage is that it is only capable of solving the minimizing compliance problem (in the current form, it cannot solve the maximizing natural frequency).

Construction of the matrix of the weight factor might be tricky. Open-source scripts usually do not support importing meshes and use a rather simple geometry. Thus, an effective script for importing meshes usually needs to be prepared. The modification of the weight factor calculation allowed us, due to the symmetry of the matrix, to solve only its upper triangle, which halved the calculating time. Another possible modification would lie in splitting the part into sub-regions, which could be calculated in parallel. Of course, the latter approach would come with its own limitations; overlapping would be necessary to be able to combine the individual regions back into the full structure, and, therefore, therefore, already solved regions would have to be calculated again.

The design performance has many criteria such as stiffness, overall solving duration, post-processing, manufacturing, etc. Each formulation has its effect on the design performance. It is difficult to pinpoint settings that are optimal for the design performance in each formulation. However, it is easy to recommend, which settings and parameter values should be avoided. It should be noted that the steps following the optimization (in particular, smoothing) might heavily affect the design performance.

It should be also noted that complex problems might need their own sensitivity analysis. However, this paper should help with the initial estimates. After finding the optimal formulations and their parameters, one could prepare scripts for automated designing of customized structures such as prosthetic aids [14,44] or scripts for automated designing of mountain climbing equipment [45].

#### **5. Conclusions**

This paper focuses on the sensitivity analysis of key topology optimization formulations. The novelty of this research comes from the presented results, which might be used in the preparation of custom scripts solving the topology optimization. The solutions were tested on various meshes with various types of elements. The paper contains important theoretical background for the problem of minimizing compliance. To have freedom in choosing such formulations, the authors prepared a MATLAB procedure solving such optimization. The prepared program allows users to import the mesh and boundary conditions. Scripts constructed within this study provided results comparable to the open-source top3d script [9]. The presented paper also includes recommendations on how to choose the parameters of topology optimization.

It is clear that uniform meshes perform generally better in this formulation; this was particularly true during optimization as it allowed the application of multiple filters.

Radius is an important part of the filtering method and correct results depend on the appropriate selection of its value. Too small a radius could possibly lead to difficult manufacturing (even if using additive manufacturing). Using a large radius could produce non-optimized shapes with grey areas. If the mesh is being refined during optimization, it is recommended to use the same or similar radius as in the previous step (coarser mesh).

The combination of the Density and Greyscale filters performed better than the Density filter alone as it yielded similar or even identical values of the objective function as the Density filter in fewer iterations. The combination also performed well on non-uniform meshes. It is obvious that the combined filter still left some grey areas but the success of this filtering was still the best of all tested filters. If the design variable in these grey areas is *x<sup>e</sup>* < 0.3, they can be removed when using penalization of *p* ≥ 3 since their contribution to stiffness is negligible.

The use of a Gaussian weight factor did not bring any advantage over the linear function as the results calculated using both functions were very similar. As constructing the matrix for the Gaussian weight factor is more difficult, a linear function is more suitable for these purposes.

In this paper, the authors used in most cases linear elements, which led to several conclusions: (i) The usage of linear tetrahedral elements is not recommended in any case. They are too stiff due to the locking effects, which greatly affects the value of the objective function. The only advantage lies in the fast calculation of displacement because the tetrahedral mesh usually has fewer nodes. However, a fine mesh would be needed to get reasonable results, which negates this only advantage. (ii) The uniform mesh provides acceptable results even if the linear approximation is used. For that reason, it is recommended to use the "Cartesian mesher" which provides uniform meshes even for complex geometries. A small disadvantage is represented by the differences in the shape of the resulting structure depending on the method of mesh creation. (iii) Quadratic elements might be less stiff but solving them would be time-consuming. If one would like to use only quadratic elements, it would be recommended to spend time preparing a better solver of linear equation systems (for example, a parallel solver).

The authors recommend performing a sensitivity analysis of the key formulations presented in the paper for each problem, regardless of whether or not the designer has previous experience with similar problems. Without the suggested analysis, more doubts arise and the creation of "false-optimal" shapes is not prevented.

**Author Contributions:** Conceptualization, M.S. and P.M.; methodology, M.S.; software, M.S.; validation, M.S., D.R. and P.M.; formal analysis, M.S. and D.R. ; investigation, M.S.; resources, P.M.; data curation, M.S.; writing—original draft preparation, M.S.; writing—review and editing, M.S., P.M., M.F. and D.G.; visualization, M.S.; supervision, P.M.; project administration, P.M.; funding acquisition, P.M., M.F. and D.G. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by The Ministry of Education, Youth and Sports from the Specific Research Project SP2021/66, by The Technology Agency of the Czech Republic in the frame of the project TN01000024 National Competence Center-Cybernetics and Artificial Intelligence and by Structural Funds of the European Union within the project Innovative and additive manufacturing technology—new technological solutions for 3D printing of metals and composite materials, reg. no. CZ.02.1.01/0.0/0.0/17\_049/0008407 and by the European Regional Development Fund under Grant No.CZ.02.1.01/0.0/0.0/15\_003/0000493 (Centre of Excellence for Nonlinear Dynamic Behaviour of Advanced Materials in Engineering) with institutional support RVO:61388998.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Data sharing not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Abbreviations**

The following abbreviations are used in this manuscript:


#### **References**


MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Symmetry* Editorial Office E-mail: symmetry@mdpi.com www.mdpi.com/journal/symmetry

MDPI St. Alban-Anlage 66 4052 Basel Switzerland

Tel: +41 61 683 77 34 Fax: +41 61 302 89 18

www.mdpi.com ISBN 978-3-0365-3176-2