MOTEO: A Novel Multi-Objective Thermal Exchange Optimization Algorithm for Optimal Design of Truss Structures

In the present paper, a physics-inspired metaheuristic algorithm is presented to solve multi-objective optimization problems. The algorithm is developed based on the concept of Newtonian cooling law that recently has been employed by the Thermal Exchange Optimization (TEO) algorithm to efficiently solve single-objective optimization problems. The performance of the multi-objective version of TEO (MOTEO) is examined through bi- and tri-objective mathematical problems as well as bi-objective structural design examples. According to the comparisons between the MOTEO and several well-known algorithms, the proposed algorithm can provide quality Pareto fronts with appropriate accuracy, uniformity and coverage for multi-objective problems.


Introduction
In the design of structures, engineers often encounter some conflicting objectives that negatively influence each other. For instance, on the one hand, they want to increase the strength of the structure to enhance safety, and on the other hand, they may want to decrease the weight of the structure for economic purposes. In mathematics, the optimization of such problems that involve two or more objectives is called multi-objective optimization problems (MOPs) [1].
Generally, there are no easy, pat answers to MOPs. Their answer is usually represented through a set of optimal solutions with a tradeoff between the optimized competitive objectives so that none of these objectives can be improved without degrading some of the other objectives. Satisfying all constraints, Pareto front provides the designer with a set of solutions. One of which can be chosen by the designer as the final design. The application of metaheuristic optimization algorithms recently became a popular research field to resolve the multi-objective design of structures [2].
Metaheuristics are mainly inspired by natural mechanisms, such as evolutionary processes, physical rules, and animals' behavior [3][4][5][6][7][8]. These approximate optimization techniques are usually able to provide acceptably good solutions for complex problems in a reasonably short time. These methods, particularly population-based metaheuristics, are efficient tools that can be applied to solve MOPs [1]. In the following, a brief overview of two popular metaheuristics, namely Genetic Algorithm and Particle Swarm Optimization, is provided.
Genetic algorithm (GA) mimics the process of natural selection using the mutation, crossover and selection operators [9]. This algorithm originally was developed for single-objective problems. NSGA II is a popular multi-objective version of GA, in which a nondominated sorting mechanism is employed to label the Pareto sets starting from the first nondominated front [10]. Afterward, a crowding distance measure is attributed to each solution using a crowded-comparison operator that guides the selection process. The algorithm prefers to survive the solutions with lower domination ranks according to the elitism notion, and it selects the solution located in less-crowded regions to preserve the diversity of solutions. Particle Swarm Optimization (PSO), as another popular algorithm, has been developed according to the swarm behavior of animals such as birds and fish [4]. In this algorithm, the particles (solutions) tend to follow the best personal experience as well as the globally best solution experienced by all particles until each iteration. In the multi-objective PSO (MOPSO) [11], a repository with a certain capacity is created in which a set of nondominated solutions are collected to be used in future steps. Furthermore, the objective function search space is divided into equal hypercubes to identify less populated places. This feature not only helps algorithm keep diversity in solutions and spread them along with the entire Pareto, but also provides designers with a wide variety of solutions rather than concentrating on specific regions.
Since structural optimization usually requires considerable computational resources to be executed, researchers constantly seek new efficient algorithms to decrease the computational costs and increase the accuracy of the solutions. Moreover, according to the no free lunch (NFL) theorem, there is no general-purpose search algorithm to outperform other methods in all problems; hence there is a strong need to develop specific modern algorithms for such specific types of problems [12]. For example, recently, Kaveh and Mahdavi [13] introduced MOCBO algorithm to minimize both displacement and weight of truss structures. Their results show that the presented algorithm can be used as an efficient technique for solving MOPs. In MOCBO, the solutions are simulated as colliding bodies, where the laws of the collision between objects are utilized for updating their variables. Luh and Chueh [14] proposed the constrained Multi-Objective Immune Algorithm (CMOIA) that considers the invading mechanism of a body against foreign substances as an analogy for the optimization process. In their research study, the cytokines concept was utilized for handling the design constraints, and the design of truss structures with displacement and weight objectives were examined. Their results indicate the effective performance of CMOIA compared to some other optimization methods. Tejani et al. [15] developed a multi-objective version of the Modified Adaptive Symbiotic Organisms Search (MOMASOS) algorithm which was developed based on the biological growth and survival process in organisms. They applied MOMASOS in structural optimization problems with mass and deflection objective functions and observed that its adaptive version provides superior results compared to those of other algorithms. In addition to the development of the new algorithms, the application of the existing algorithms to various MOPs also engrossed the interest of many researchers. For example, Nan et al. [16] recently optimized the coordinate of the node space truss structures using a multi-objective evolutionary algorithm (MOEA). Any design with suitable boundary conditions and topology can be selected out of their final Pareto set providing designs with both optimum structural weight and minimum displacement. Other similar researches on the development and application of multi-objective optimization can be found in Refs. [17][18][19][20][21][22].
In the present paper, a multi-objective version of the newly-established physics-based algorithm, socalled TEO [23,24], is introduced, and it is applied to several structural examples with weight and displacement objective functions. Various mathematical benchmark functions are also examined to show the accuracy and quality of its Pareto front compared to true Pareto. The results are compared with those of some well-known algorithms.
The remainder of the paper has been arranged as: In Section 2, the basic concepts of multi-objective optimization are presented; in Section 3, the Thermal Exchange Optimization and its multi-objective version (MOTEO) is presented; in Section 4, the results of mathematical and structural examples are discussed; and the concluding remarks are finally provided in Section 5.

Basic concepts in multi-objective optimization
This section encompasses some necessary definitions for multi-objective optimization theory that will be used in the next sections. In these definitions, the optimization problems are assumed to be minimization problems.

Definition 1. The multi-objective optimization problem (MOP)
A MOP is defined as follows where, ( ) denotes a set of at least two cost functions ( ( )), which should be minimized; is the space of feasible solutions = ( 1 , 2 , … , ), which are usually constrained between maximum and minimum bounds corresponding to each variable. MOPs often possess more than one solution, none of which can dominate the others. The dominance concept is defined in the next paragraph.

Thermal Exchange Optimization
Physics-inspired algorithms employ physical phenomena to update search agents inside the search space and thereby solving optimization problems. For example, Gravitational Search Algorithm (GSA) uses Newtonian gravity law, motion laws, and mass interactions to move the search agents (masses) [25]; and Big Bang and Big Crunch (BB-BC) algorithm mimics a theory with the same name to distribute the search particles randomly throughout the search space in the Big Bang phase and shrinks them to a single point as the center of mass [26]. Recently, Kaveh and Dadras [24] proposed a new metaheuristic optimization algorithm, so-called Thermal Exchange Optimization (TEO), which was inspired by Newtonian cooling law. In this algorithm, the search agents are divided into two searching groups, including environment and cooling objects whose temperature represents the optimization variables. In the following, firstly the steps of the original TEO algorithm for a single-objective problem are stated, and the necessary changes for devising a multi-objective version are then explained.

Step 1. Initialization of temperatures
First, the search agents-temperatures of objects-are randomly initialized throughout the search space. These solutions can be produced by Eq. (5).
where is the number of search objects; denotes the initial solution of the th object; and are the vector of upper and lower bounds of solutions, respectively; and is a vector of random numbers, independently generated for th object, with a uniform distribution between 0 and 1.

Step 2. Evaluation of objects
After assigning temperatures to the objects, they are evaluated by cost functions. The objects are sorted according to their cost function values (costs) in descending order, and the first number of objects-equal to the assumed number of objects-are preserved.

Step 3. Saving the best objects
Some of the historically best solutions are saved in Thermal Memory (TM). Therefore, in each iteration, if the new better solution(s) with a lower cost than the current TM is found, TM is updated and the new solutions with lower costs are saved. This step helps the algorithm keep the position of promising regions in the search space and boost the exploitation (intensification).

Step 4. Creating groups
The sorted objects are paired, such that the first half is considered as the environment for objects in the second half and vice versa. For example, as shown in Fig. 2, 1 is paired with 2 +1 .

Steps 5 and 6. Calculating and parameters
Newton's law of cooling states that the heat loss rate of a lumped object is directly proportional to the difference in the temperatures between the body and its surroundings, as formulated in Eq. (6).
where ( ) is the new temperature in time after the thermal exchange between an object with a temperature 0 and the environment with temperature , and is a constant relating to several parameters such as the object's specific heat capacity and density. As can be seen, when has a higher value, the object tends to change less amount of temperature. In TEO, inspired by this phenomenon, is defined such that the less the cost, the less the changes in the solution. Therefore, is defined as follows = where, and respectively are the cost of the current object and the worst object with the highest cost among the population. Also, the time parameter is related to the iterations and is defined as Eq.
where denotes current iteration number and stands for the maximum permitted iterations.
Step 7. Escaping local optima (i) In order to escape from local minimums, the metaheuristics should regulate the solutions employing some stochastic mechanisms. For instance, GA uses mutation to avoid getting trapped in such solutions and to explore various chromosomes. In TEO, the current step is one of the twin mechanisms that are dedicated to this task. Here, the environmental solutions are randomly modified before temperature updating using Eq. (9).
and are the temperatures of the object before and after the modification; 1 and 2 are internal controlling parameters; and is a random vector containing random numbers between 0 and 1. The equation has been devised so that it can reduce the randomness as the algorithms reach the last iterations and thereby reducing the exploration and increasing the exploitation.
Step 8. Temperature updating According to Eq. (6), the new temperature of cooling objects can be calculated as follows.
and are the new updated and previous temperature of the th object, and the and parameters were defined in Steps 5 and 6.

Step 7. Escaping local optima (ii)
This step is the second mechanism designed to avoid trapping in local minima so that TEO can explore the global optimum solution. Herein, with the probability of one of the dimensions is randomly selected and its value is changed to a random number between its allowable upper and lower bounds.

Step 8. Termination conditions
The maximum number of iterations is controlled in this step, and if the condition is met, the algorithms cease the search process and report the best solution found so far; otherwise, it resumes and returns to Step 2.

Multi-objective thermal exchange optimization
The steps of the multi-objective Thermal Exchange Optimization (MOTEO) version are provided in this section. Since MOTEO has structural similarities with the single-objective version, for sake of brevity just different steps are mentioned. The differences lie in the way that the objects are sorted and the way that parameter is measured.

Selection mechanism
In multi-objective problems, dependent on the number of objectives functions, several cost values are attached to each object. Therefore, an idea regarding their ranking is the sorting based on their Pareto fronts. In this regard, a non-dominated sorting algorithm is executed to divide the solutions into several fronts [10]. A non-dominated sorting algorithm can be designed as follows: an index (equal to one) is assigned to each object, and the cost values are compared with those of other objects, as stated in Section 2. Any time that the object is dominated by another object, its index value is increased by one. After performing this step for all objects, the objects with the same index values are considered as a separate group. It is apparent that the objects' index values have an inverse relationship with their quality.
After performing the non-dominated sorting algorithm, the rank of each group against other groups is determined. However, their internal ranking has not been specified yet. In order to compare the objects inside a group, their crowding distance, as a secondary criterion, is calculated and compared with other groups. Crowding distance is calculated by Eq. (11) [10] where, is the crowding distance of the th object; denotes the distance of two neighboring solutions of the th object considering the th cost function. This parameter for a bi-objective problem is illustrated in Fig. 3. and are respectively the cost of the best and worst objects considering the th cost function. It is noteworthy that the crowding distance of the best and worst objects is assumed as infinity.

Fig. 3. An illustration of distance parameter
To preserve diversity, the crowding distance of the object located inside each group are compared, and those having higher distance attain superior ranks. To put it simply, if a solution belongs to a better front than , the former receives a better rank than the latter. If the front ranks are the same but > , receives a relatively better rank because of higher crowding distance.
The selection mechanism utilized in MOTEO is pictorially illustrated in Fig. 4. As seen, in each iteration the number of populations will remain constant.

parameter
According to the ranking system explained above, is redefined as follows: where, is the final rank of the solution according to its Pareto front and crowding distance, and stands for the number of populations.
The rest steps of the MOTEO are the same as TEO, and MOTEO can be implemented only by making the mentioned adjustments, as demonstrated in Fig. 5.

Numerical Examples
Eight mathematical test functions and two well-known truss structures will analyze the efficiency of the proposed method. The structural problems and mathematical functions are employed to verify the ability of the multi-objective optimizers to handle problems with various characteristics such as non-convexity and non-linearity. The algorithm has been coded in MATLAB 2020b and the structures are analyzed with the direct stiffness method. The present work is performed with the following features on the computer: CPU 2.3 GHZ (an Intel Core i9 computer platform), Ram 16 GB 2400 MHz DDR4 on a computer with Macintosh (macOS BigSur).

Performance Metrics
We used four metrics in the following way to assess the results of the algorithms. [27,28]:

Generational distance (GD):
This index calculates the distance between the true Pareto ( ) and the obtained Pareto front ( ). This distance is calculated by Eq. (13).
where is the number of non-dominated solutions in and is the Euclidean distance between the solution in and its nearest solution in . It should be noted that a smaller value of implies a more desirable convergence of the algorithm.

Spacing (S):
This measure indicates how evenly the obtained solutions are distributed along the . This measure is formulated as follows: where is the number of non-dominated solutions in and is the Euclidean distance between the solution in and the nearest solution in . The smaller value of shows better convergence.

Maximum Spread (MS): This metric evaluates how ''well'' the
is covered by the through hyper-cubes formed by the extreme cost values observed in the and . It is defined as: where is the number of objectives, and are the maximum and minimum of the ℎ objective in , respectively; and and are the maximum and minimum of the ℎ objective in , respectively. A larger value of MS indicates a better spread of solution as well.

Inverted Generational Distance (IGD):
It is a metric for assessing the quality of approximations to the Pareto front obtained by multi-objective optimization algorithms. This measure is defined as: here, is the number of elements in , and is the Euclidean distance between the solution in and its nearest solution in . The Zero value shows that all of the generated elements are in the true front of Pareto. [29]

Mathematical Test Function
In this section, the results of the proposed multi-objective optimization algorithm have been analyzed with eight well-known benchmark problems. [27,28]. These problems are introduced in the following and their formulations are provided in Table 1.

Table 1
Multimodal benchmark functions with fixed-dimension Function d Range  The results of the presented method are compared with the three famous multi-objective optimization algorithms-Strength Pareto Evolutionary Algorithm two (SPEA2) [28], Non-dominated Sorting Genetic Algorithm II (NSGA-II) [10], and Multi-Objective Particle Swarm Optimization (MOPSO) [11] algorithms. The internal parameters considered for these algorithms are shown in Table 2. In this study, the number of function evaluations for all algorithms is set to 300000 and the results are reported based on 30 independent runs.
The mean and standard deviation of GD, S, MS, and IGD metrics for bi-objective test functions along with tri-objective problems are presented in Tables 3-10. It should be mentioned that the boldfaced values indicate the best results obtained among the compared algorithms. As seen, MOTEO has the best performance in terms of GD, S, MS and IGD metrics. Therefore, it can be said that, compared to other algorithms, MOTEO has the shortest distance from true Pareto fronts, provides a uniform distribution of solutions, and covers a wider range of cost values. In terms of GD, S, and IGD, the NSGA-II is only second to MOTEO. NSGA II. The obtained best Pareto fronts of MOTEO for five bi-objective and three tri-objective problems are presented in Figs. 6 and 7, respectively. The figures prove that MOTEO is able to find a Pareto optimal set near the global Pareto front and to preserve solution diversities across space, as well. All in all, it can be inferred that most results of MOTEO are superior to those of other algorithms.

Table 5
The mean (and standard deviation) results of benchmark functions MS performance metric

Truss Structure Examples
In where { } is the vector of design variables, is the number of design variables, the structure weight is described as ({ }), the displacement is described as ({ }), and the number of constraints is denoted as . Here, the th member's material density is shown by , is its length, is the member's cross-sectional area, and ({ }) denotes design constraints. Also, min and max are the lower and upper bounds of the design variable vector, respectively. As the penalty function is simpler and easier to carry out, the limitations expressed in the following are handled: where is the sum of the violation ratios of constraints and constants 1 and 2 are selected to enhance the exploration and exploitation rate of the search. Here, 1 is set to 1, and 2 is set to 1.5 at the first iteration of the search process and is gradually increased to 3 [30].
The allowable compressive and tensile stress constraints according to the AISC [31] is considered as: where the modulus of elasticity is described by , the slenderness ratio is shown by , is defined as √2 2 / , and the yield stress of steel is .

A 120-bar Dome Truss
The optimization of 120-bar dome truss has been studied by different researchers to minimize the weight as a single objective [32,33]. In this paper, the total structural weight as well as the maximum displacement in all directions on all the nodes are considered as two optimization objectives. The density of the material and the modulus of elasticity are respectively 0.288 / 3 (7971.810 / 3 ) and 30.450 (2.1 × 10 5 ), and the yield stress is taken as 58 (400 ). This structure is drawn in Fig. 8. The dome is subjected to the vertical loading through all free nodes as following: −13.49 (−60 ) at node 1, −6.744 (−30 ) at nodes 2 to 14, and −2.248 (−10 ) at the remaining nodes. The truss members are categorized into seven independent size variables, and the cross-sectional areas are between 0.775 2 (5 2 ) and 20 2 (129.032 2 ).
As the true front of Pareto is not correctly identified, here the performance measurement used in the previous section cannot be applied. However, the best two extreme-point values in Table 7 are compared in order to provide a quantitative assessment of the performance of the optimizers. The extent of the obtained non-dominated front should be maximized, i.e., for each objective, a wide range of values should be covered by the non-dominated solutions [28]. As provided in the table, MOTEO is able to find solutions with lower displacement and weight in the extreme points. The best Pareto fronts of NSGA-II, SPEA2, MOPSO, and MOTEO obtained after 20 runs are plotted and compared in Fig. 9. It can be seen that the obtained set of solutions for MOTEO is distributed uniformly. Additionally, the obtained Pareto front by MOTEO generally dominated those obtained by MOPSO, SPEA2 and NSGA-II. The majority of the Pareto front can also be covered in a suitable convergence.

A 582-bar Tower Truss
The single-and multi-objective optimization of the 582-bar tower truss structure, shown in Fig 10, has been studied in the literature [34]. Here, this structure with continuous sizing variables as a multiobjective optimization problem is investigated. Two objectives, namely the total used steel volume and the maximum displacement in all degrees of freedom on all nodes are optimized simultaneously. The members are symmetrically divided into 32 groups. The truss is subjected to the lateral loads in both xand y-directions and a vertical load in the z-direction at all free nodes as following: 1.12 (5 ), 1.12 (5 ) and −6.744 (−30 ), respectively. The lower and upper bounds of cross-sectional areas are taken as 3.1 2 (20 2 ) and 155 2 (1000 2 ). The yield stress of steel is 36 (253.1 ) with the elasticity modulus of 29000 (203893.6 ). According to ASD-AISC, the slenderness ratio for tension and compression members are respectively limited to 300 and 200 [31].
(a) (b) Fig. 10. Scheme of the 582-bar tower truss structure: (a) perspective view (b) dimensions. Fig. 11 depicts the Pareto fronts obtained by NSGA-II, MOPSO, SPEA2, and MOTEO algorithms. As seen, the Pareto front obtained by MOTEO generally dominates the rest of the algorithms. Like the preceding example, the true Pareto front is not known for this example, and the performance metrics are not applicable. Table 8 presents the two best extreme point values obtained by the algorithms proving that MOTEO covered a larger portion of the Pareto front. As can be seen, the minimum of both displacement and volume is obtained by MOTEO. These results show the proper convergence and sparsity efficiency of the MOTEO in optimizing problems with multiple conflicting functions.

Conclusion
In this paper, a multi-objective version of the recently-developed Thermal Exchange Optimization is proposed. The differences between the introduced Multi-objective Thermal Exchange Optimization (MOTEO) and the basic single-objective version mainly rely upon the definition of parameter and ranking scheme. The rank of solutions is determined using a non-dominated sorting mechanism and crowded distance between obtained solutions, as illustrated in Fig. 4. According to the new ranking system, parameter is calculated for each solution, and the position of solutions is updated accordingly.
In order to investigate the proposed MOTEO, firstly eight bi-and tri-objective mathematical problems have been selected and MOTEO's performances are compared to three well-established algorithms, in terms of four accuracy and Pareto quality metrics. The results reveal the efficiency and competitive performance of the MOTEO considering most of the measures and criteria. Also, to examine MOTEO in the optimization of real-world engineering problems, two structural examples with higher dimensions have been taken into account. It was observed that MOTEO is capable of solving multi-objective structural problems effectively and demonstrates a competitive and satisfactory performance compared to other well-known algorithms.