Optimal Design of Electrothermal Microactuators for Many Criteria by Means of an Immune Game Theory Multiobjective Algorithm

: The paper presents the application of the IMGAMO (immune game theory multiobective algorithm) in the optimal design of electrothermal microactuators. Several numerical tests on the mathematical benchmark test functions were performed, showing the superiority of the IMGAMO, when a large number of criteria are considered, in comparison to other multiobjective optimizers. A parametric numerical model of an electrothermal microactuaror was developed and veriﬁed. Six functionals, which depend on various thermal and mechanical quantities of the microactuator, were proposed, formulated and numerically implemented. These functionals represent real requirements asked of microactuators. The boundary-value problem of an electro-thermo-mechanical ﬁeld was solved multiple times during the course of optimization by way of the ﬁnite element method (FEM). A numerical example of multiobjective optimization of chevron-type electrothermal actuators is included in the paper. Representation of the multi-dimensional Pareto fronts by means of scatter plot matrices, aided by self-organizing maps (SOMs), is presented. The novel method of selecting interesting, compromise-solutions is proposed and described.


Introduction
Electrothermal microactuators are mechanical systems where thermally-induced expansion and contraction of materials is responsible for the induction of motion. A temperature change inside the actuator is caused by Joule heating due to the electrical current flow. In the context of microscale and nanoscale engineering, thermal actuators can be used as microelectromechanical systems (MEMS) capable of translating electrical signals into controllable displacements and forces [1]. Electrothermal microactuators have proven to be very useful MEMS structures, which are commonly used in telecommunication, automotive, aviation and medical industries; e.g., microgrippers, switches, filters and modulators. Compared to the other types (e.g., magnetostatic, piezoelectric and electrostatic), the electrothermal actuators have many advantages. They can generate large forces or displacements per unit volume; require low actuation voltages; have linear material characteristics (strain-stress relationship), despite the high temperatures occurring in them; and moreover, they are relatively easy and cheap to fabricate. The disadvantages of microthermal actuators concern comparatively high-power consumptions and relatively slow response times. The Joule heating effect is possible due to the high value of electrical resistance of the material, which is the polycrystalline silicone. Mainly, two types of such microactuators are used in MEMS: U-beam and V-beam ones. U-beam type actuators were originally demonstrated in the early 1990s, whereas V-beam ones were in the late 1990s, and quickly became very popular for the generation of mechanical motion. U-beam microactuators consist of two beams (arms) which have different cross-section areas. They produce different amounts of Joule heating in the presence of the electric currents, which causes the bending of the apexes of the actuators. For the V-beam microactuator, arms are fabricated with a pre-bending angle. The longitudinally extending arms, also due to the Joule heating effect, move the central shaft. The central shaft aggregates the actuation force of the multiple beams, resulting in a total actuation force [2]. The V-beam actuator, called variously-chevron type, bent-beam or symmetric thermal actuator, is considered in the paper. Figure 1a presents an example shape of the chevron microactuator consisting of four pairs of arms, whereas Figure 1b shows the deformation and map of equivalent stress in the structure.
Appl. Sci. 2019, 9, x FOR PEER REVIEW 2 of 29 actuators were originally demonstrated in the early 1990s, whereas V-beam ones were in the late 1990s, and quickly became very popular for the generation of mechanical motion. U-beam microactuators consist of two beams (arms) which have different cross-section areas. They produce different amounts of Joule heating in the presence of the electric currents, which causes the bending of the apexes of the actuators. For the V-beam microactuator, arms are fabricated with a pre-bending angle. The longitudinally extending arms, also due to the Joule heating effect, move the central shaft.
The central shaft aggregates the actuation force of the multiple beams, resulting in a total actuation force [2]. The V-beam actuator, called variously-chevron type, bent-beam or symmetric thermal actuator, is considered in the paper. Figure 1a presents an example shape of the chevron microactuator consisting of four pairs of arms, whereas Figure 1b shows the deformation and map of equivalent stress in the structure.
(a) (b) The process of design of such structures has to incorporate many technical aspects and restrictions. Accordingly, optimal design of electrothermal microactuators, requires consideration of many criteria, which depend on the quantities like actuation force, generated displacement, heat or electrical losses, stress concentration, risk of buckling and dynamic characteristics. For such structures, the above criteria strongly depend on the geometry. These criteria are very often contradictory, which leads to the necessity of accepting compromise-solutions. Analytical formulas related to basic electrical, thermal and mechanical quantities have been derived only for the simplified geometry (straight shape arms) of the microactuators. For the arbitrary shape of the geometry in order to calculate objective functions, the quantities derived from numerical simulations of coupled electro-thermo-mechanical analysis have to be used. Functionals expressing assumed criteria, calculated on the basis of multiphysics numerical simulations can be strongly multimodal (have many local minima). It significantly limits the scope of the applicability of classical gradient or non-gradient optimization methods. It can be seen that in order to efficiently optimize such structures, proper multiobjective optimization techniques have to be used [3]. Application of the multiobjective optimization methods leads to a situation where the result is not only a single solution, but a number of optimal solutions. These solutions are called Pareto-optimal solutions.
Very popular approaches to multiobjective optimization are the weighted sum method or taking only one criterion into account, while imposing other criteria as constraints. Such methods are handy to use, because they do not require modification of the optimization algorithm (the single-objective optimization algorithm can be used). Unfortunately, they have significant limitations, which reduce the scope of their application, such as: • Difficulties in selection of the balance between criteria before optimization, having no knowledge of the considered problem; • Certain assumptions imposed before optimization can lead to improper or unphysical results; The process of design of such structures has to incorporate many technical aspects and restrictions. Accordingly, optimal design of electrothermal microactuators, requires consideration of many criteria, which depend on the quantities like actuation force, generated displacement, heat or electrical losses, stress concentration, risk of buckling and dynamic characteristics. For such structures, the above criteria strongly depend on the geometry. These criteria are very often contradictory, which leads to the necessity of accepting compromise-solutions. Analytical formulas related to basic electrical, thermal and mechanical quantities have been derived only for the simplified geometry (straight shape arms) of the microactuators. For the arbitrary shape of the geometry in order to calculate objective functions, the quantities derived from numerical simulations of coupled electro-thermo-mechanical analysis have to be used. Functionals expressing assumed criteria, calculated on the basis of multiphysics numerical simulations can be strongly multimodal (have many local minima). It significantly limits the scope of the applicability of classical gradient or non-gradient optimization methods. It can be seen that in order to efficiently optimize such structures, proper multiobjective optimization techniques have to be used [3]. Application of the multiobjective optimization methods leads to a situation where the result is not only a single solution, but a number of optimal solutions. These solutions are called Pareto-optimal solutions.
Very popular approaches to multiobjective optimization are the weighted sum method or taking only one criterion into account, while imposing other criteria as constraints. Such methods are handy to use, because they do not require modification of the optimization algorithm (the single-objective optimization algorithm can be used). Unfortunately, they have significant limitations, which reduce the scope of their application, such as:

•
Difficulties in selection of the balance between criteria before optimization, having no knowledge of the considered problem; • Certain assumptions imposed before optimization can lead to improper or unphysical results; • Performing a single optimization task gives only one solution; • The impossibility to obtain solutions for non-convex Pareto fronts.
Population-based approaches, such as metaheuristic algorithms, are free from the disadvantages mentioned above. Moreover, because such algorithms process not one, but a set of solutions, and the result of multicriteria optimization is a set of optimal solutions in the sense of Pareto, they are an ideal tool for these type of tasks.
A group of papers are devoted to the shape optimization of the different types of the thermal actuators: Lee et al., Huang et al. and Ionescu perform optimization of the characteristic dimension of the laterally driven microactuator [4][5][6]. Kwan et al. proposed optimization of the chevron type microactuator, assuming nonuniform cross section of the beams [7]. The process of optimal design can be realized assuming the variant analysis [8,9] or by means of different optimization methods [10]. Some papers are devoted to the topological optimization, which allows finding desired shape of the actuator with respect to applied restriction and criteria. Usually, the optimal shape of the actuator consists of a set of variously shaped arms in a number of geometrical configurations. Sigmund, Giusti et al. and Ansola et al. demonstrated maximizing the output displacement in a given direction of the actuator by means of topology optimization, obtaining quite complicated optimal shapes of the actuator [11][12][13]. Furuta et al. and Xia et al. demonstrated the topological optimization for thermoelectric actuators aided by level-set methods [14,15]. The result of topological optimization can be a good tip for the designer on what the best configuration of the arms is with regard to the considered criterion; however, very complicated shapes (common for topological optimization) can be an obstacle in the design process of practical shapes consisting of repetitive pairs of arms (which is needed for aggregation of the actuation force in the central shaft by multiple arms). Heo and Kim present topological optimization for the rotary thermal actuators [16]. For the majority of cases, a process of finding an optimal design of thermal actuators is dealt with as a single optimization problem. Multiobjective optimization methods are not popular in the design of electrothermal actuators; however, some papers have been devoted to it. Di Barba et al. solved the bi-objective optimization problem of a U-shaped actuator by means of the NSGAII algorithm (non-dominated, sorting genetic algorithm) [17]. It is a well-known fact that NSGAII algorithm (probably the most popular multiobjective optimization algorithm used in diverse engineering disciplines) does not work well for more than three criteria (practically 2). For the other types of the actuators (e.g., piezoelectric and electrostatic) multiobjective optimization problems generally also deal with the simultaneous consideration of at most two criteria.
A review of the literature allows us to claim a lack of research that incorporates many criteria into the optimal design of thermal actuators. Moreover, shape optimization tasks are devoted to the actuators with simplified geometry of the arms. More sophisticated shape of the arms in the optimal design has been considered in the case of topological optimization, but the disadvantages of such an approach have been described above. Furthermore, for the topological optimization, only single objective optimization tasks have been considered.
Bi-objective optimization of the U-beam actuator has been presented in [18], whereas in the present paper, a novel and efficient method of multiobjective optimization for higher number of criteria is proposed. The proposal combines an algorithm based on artificial immune systems (AIS) and elements of game theory (GT); it has been developed and implemented by authors.
Game theory deals with the analysis of conflict situations and can be applied to problems in fields such as economics, psychology and recently, engineering as well. The work which is assumed to be the beginning of game theory is accredited to John von Neumann and Oskar Morgensern [19]. In game theory, the conflict situations are named games. They have several essential features:

•
Number of players is known and finite; • Full game rules are set; • Payoffs are determined by payoff functions, which depend on strategies chosen by players.
A strategy is a complete description of a player's activity in all situations. Ideas of multiobjective optimization (MOO) and game theory have several common features. The game theory problems were considered applying multiobjective optimization concepts. MOO solutions which best take into account all criteria (so called non-dominated solutions) are sought, and in game theory, players are looking for best strategies. Finally, in game theory we try to find equilibrium and in MOO we try to find a set of solutions optimal in a Pareto sense. Several papers which considered these approaches were published [20][21][22][23][24].
Artificial immune systems, on the other hand, are one of the systems which belong to a wide range of methods of computational intelligence. The most developed algorithms are: immune networks, clonal selection and negative selection. Clonal selection, used in this present paper, is dedicated to solving optimization problems: fitness functions correspond to the affinities of antibodies to antigens. Solutions in a population are represented by antibodies. Each antibody is cloned several times like in the real immune system. The number of clones is proportional to an antibody's affinity to antigen (its value of the fitness function) and then the hypermutation is applied. These ideas are extended with the suppression that takes place [25]. Immune systems and algorithms succeed in their application to global and multimodal optimization [26]. There are several adaptations and modifications of the clonal selection algorithm applied to problems with many criteria [27][28][29][30][31].
Artificial immune systems, compared to the population-based metaheuristic algorithms, seem to be predisposed to solve multiobjective problems. The main features of them are directly inspired by the immune system. Diversity of solutions is preserved thanks to a suppression mechanism, which eliminates similar and useless solutions, and as a result, controls the size of population. The elitism mechanism helps preserve the best solutions during the optimization course.
When more than two criteria are considered, the proposed multiobjective optimization algorithm significantly outperforms other metaheuristic multiobjective algorithms (NSGAII, SPEA2) [17,32]. In the present paper, six objective functionals have been proposed, formulated and implemented. The shape optimization method proposed in the paper allows the design of the arms of the chevron type microactuator with the smooth nonuniform cross sectional area and smooth transition between arms and anchors or central shaft. Non uniform, rational B-splines (NURBS) have been applied for the parametrization of the arms [33]. It is obvious that stress concentration is located in sharp corners (see Figure 1b). Increasing the radius between actuator beams and the anchors or central shaft can significantly reduce the stress concentration. Figure 2 shows a comparison of the distribution of the equivalent von Mises stress between a geometric arrangement with a sharp connection of the arms and one with rounded corners (transition radius is twice larger than the width of the arm) in the actuator for the same boundary conditions. This example shows that the maximal value of the equivalent stress is significantly reduced from 982 MPa to 476 MPa. The proposed method of optimal design allows for both, finding the arbitrary shape of the arms and precise selection of the radius at the rounded corners, taking into account the optimal selection of the pre-bending angle as well. As mentioned previously, the result of the multiobjective optimization is expressed in not one, but a set of optimal solutions. When a large number of criteria is considered, the process of selection of a particular compromise-solution could be problematic. The alternative method, compared to the traditional one, concerning the representation and selection-compromise-solutions, is proposed and demonstrated in the paper. The proposed method is aided by self-organizing maps (SOMs), which is one type of artificial neural network (ANN).
pre-bending angle as well. As mentioned previously, the result of the multiobjective optimization is expressed in not one, but a set of optimal solutions. When a large number of criteria is considered, the process of selection of a particular compromise-solution could be problematic. The alternative method, compared to the traditional one, concerning the representation and selection-compromisesolutions, is proposed and demonstrated in the paper. The proposed method is aided by selforganizing maps (SOMs), which is one type of artificial neural network (ANN). The rest of the paper is organized as follows. In Section 2, the model of the chevron type actuator is described, including the governing equations of the problem considered, the numerical model developed, and its verification. Section 3 contains the formulation of the multiobjective optimization problem, including the definition of the formulated criteria and SOMs, whereas in Section 4 the system of optimization build by means of IMGAMO algorithm is presented. Section 5 presents details of the optimization model of the electrothermal microacuator. On the basis of the research carried out, in Section 6, results and discussion are presented. The work is concluded in Section 7 with final remarks. Appendix A shows the details of IMGAMO algorithm, whereas Appendix B illustrates its effectiveness in comparative tests with other algorithms for several benchmarks.

Governing Equations
As mentioned in the previous paragraph actuation of the thermal microactuator is caused due to the electrical current I, which is converted to the heat according to the Joule-Lentz law [34]: where Q is internal heat generated, R is electrical resistivity and t is the time.
The thermal strain occurring in the longitudinal arms causes structural deflection. Usually such an actuator is affected by the voltage of the magnitude of a few volts, which causes an increase of temperature, even above 1200 K. Nevertheless, the linear behavior of the microactuator with respect to electric, thermal and mechanical effects is preserved. Such a phenomenon is described by the partial differential equation of electrostatics, Equation (2); heat conduction, Equation (3); and thermoelasticity, Equation (4) [35].
where φ is electric potential, T is the temperature, u represents the displacement values, ρ is charge flux density, ε 0 is electric constant, k is thermal conductivity, Q is internal heat source, α t is the linear expansion coefficient, and µ and λ are the Lamé constants, which can be expressed as follows: where ν, E and G are Poisson's ratio, Young's modulus and the shear modulus, respectively.
The structure considered here is an electro-thermo-elastic body of the domain Ω, with the boundary area Γ (Figure 3). Partial differential Equations (2)-(4) have to be supplemented by boundary conditions. where and are known temperature and heat flux on the parts T and q of the boundary. The third type of thermal boundary condition (Robin condition) is convection condition, for which ( ) T x ∞ is an ambient temperature, whereas α is a heat convection coefficient.

•
For the elasticity problem boundary condition can be defined as: where 0 ( ) u x and 0 ( ) p x are known displacement and mechanical loads on the parts u Γ and p Γ . Separate parts of the boundaries where electrical, thermal and mechanical boundary conditions are specified define the boundary of the body. When no electrical, thermal or mechanical boundary conditions are specified on the particular parts of the boundary, it corresponds to the charge flux density free, heat flux free or traction-free conditions respectively. The relationships between separate parts of the boundary Γ can be written as: The finite element method (FEM) has been used for solving the electro-thermo-mechanical problem [36,37]. The partial differential equation after discretization and after taking into account boundary conditions, can be written as a set of algebraic equations and given in the following matrix form: • For electrostatics: where φ 0 (x) and ρ 0 (x) denote known electric potential and electric charge flux density on the appropriate parts of the boundary Γ φ and Γ ρ ; • For the heat conduction problem, boundary conditions are as follows: x ∈ Γ T : where T 0 (x) and q 0 (x) are known temperature and heat flux on the parts Γ T and Γ q of the boundary. The third type of thermal boundary condition (Robin condition) is convection condition, for which T ∞ (x) is an ambient temperature, whereas α is a heat convection coefficient.
• For the elasticity problem boundary condition can be defined as: where u 0 (x) and p 0 (x) are known displacement and mechanical loads on the parts Γ u and Γ p . Separate parts of the boundaries where electrical, thermal and mechanical boundary conditions are specified define the boundary of the body. When no electrical, thermal or mechanical boundary conditions are specified on the particular parts of the boundary, it corresponds to the charge flux density free, heat flux free or traction-free conditions respectively. The relationships between separate parts of the boundary Γ can be written as: The finite element method (FEM) has been used for solving the electro-thermo-mechanical problem [36,37]. The partial differential equation after discretization and after taking into account boundary conditions, can be written as a set of algebraic equations and given in the following matrix form: The global electrical conductivity matrix K E , global thermal conductivity matrix K T and global stiffness matrix K M are assembled over element matrices. V,I,T,Q,U and F are the nodal vector of voltage, electric current, temperatures, heat fluxes, displacements and applied forces, respectively. The problem is weakly coupled, so electrical, thermal and mechanical analyses can be solved separately. Coupling is carried out by transferring loads between the analyses considered and by using staggered procedures. After solving Equation (10), the nodal vector of heat generation due to the current flow Q E is calculated, whereas after solving Equation (11), the nodal vector of forces due to the thermal strain vector F T is calculated.
For the buckling analysis, the finite element formulation takes the form: Linear buckling analysis is assumed. Equation (13) describes the eigenvalue problem, which is solved by the Lanczos method. Matrix ∆K G is a function of the load increment ∆P; λ i is the eigenvalue, whereas eigenvector v i allows for the calculation of post-buckling deformation mode.
As mentioned in the first paragraph, one of the criteria formulated depends on the contact force. Because the movement of the central shaft is almost perfectly vertical, friction between the end of the central shaft and the reaction plane can be neglected. The contact problem is modelled as rigid, deformable and frictionless. The thermal actuator is modelled as a deformable body, whereas the rigid surface is placed with the arbitrarily chosen initial gap between the tip of the central shaft and the rigid surface. For the contact between a deformable body and the rigid surface, the constraint responsible for maintaining the lack of penetration is implemented by transforming the degrees of freedom of the contact node and applying a boundary condition to the normal displacement.

Verification of the Numerical Model
The optimal design of the chevron type thermal actuator requires building a numerical model of the actuator, which allows one to calculate the values of criteria with the satisfactory level of accuracy. It is possible to analytically calculate electrical, thermal and structural quantities for the actuator with straight arms. This subsection is dedicated to the verification of the numerical model prepared. Simplifying the verification procedure comparison between the numerical and analytical models only for mechanical quantities is shown. Since the thickness perpendicular to the midplane of the actuator (z-axis) and boundary conditions are constant, and thickness is very small compared to other dimensions in the model, only two-dimensional models will be considered.
Consider the fundamental unit of the V-type actuator, a pair of angled beams ( Figure 4). The average actuator temperature rise ∆T, caused by the flow of electric current, produces the output displacement u y and the force F y . An external force F i may be applied to the central shaft in y direction, affecting the u y displacement. Analytical modeling of the V-type actuators can be found in [2,5]. The equation below shows relationships between those quantities.
If the force constraint F i is equal to zero, Equation (14) expresses the relationship for the free-moving actuator displacement. Taking into account-average temperature rise ∆T = 450 K, E = 158, 000 MPa, length of the beam L = 220 µm and the square cross-section area A = 25 µm 2 , the FEM numerical model in MSC Mentat/Marc software was prepared and compared with the analytical solution. Plane stress model consisting of four node quad elements was assumed. The density of the mesh was established after some preliminary calculations. The number of 6-8 elements per width of the beam ensures correct results when the actuator bends, and therefore, the model consists of about 6000 elements ( Figure 5).
Appl. Sci. 2019, 9, x FOR PEER REVIEW 8 of 29 the beam ensures correct results when the actuator bends, and therefore, the model consists of about 6000 elements ( Figure 5).   Figure 6 shows the characteristics of the output displacement with respect to the pre-bending angle. The numerical simulation results are in close agreement with the analytical results. The highest discrepancy was obtained for the peak value of the displacement, where relative error was equal to 2.3%.    Figure 6 shows the characteristics of the output displacement with respect to the pre-bending angle. The numerical simulation results are in close agreement with the analytical results. The highest discrepancy was obtained for the peak value of the displacement, where relative error was equal to 2.3%.  Figure 6 shows the characteristics of the output displacement with respect to the pre-bending angle. The numerical simulation results are in close agreement with the analytical results. The highest discrepancy was obtained for the peak value of the displacement, where relative error was equal to 2.3%.
Appl. Sci. 2019, 9, x FOR PEER REVIEW 8 of 29 the beam ensures correct results when the actuator bends, and therefore, the model consists of about 6000 elements ( Figure 5).   Figure 6 shows the characteristics of the output displacement with respect to the pre-bending angle. The numerical simulation results are in close agreement with the analytical results. The highest discrepancy was obtained for the peak value of the displacement, where relative error was equal to 2.3%. Analytical Equation (14) may be useful for the straight arms and basic design only; moreover, it underrepresents the interaction between constraint force and displacement when the displacement is close to zero or when a large force is applied [2]. From the practical point of view, actuation force should be modeled as a contact force, as in the model in [38]. Figure 7 shows the relationship between the pre-bending angle and the contact force calculated at the rigid surface, for the same model of the actuator. The gap is understood as a distance from the tip of the central shaft to the rigid surface. A small value of the pre-bending angle actuator cannot generate a high enough value of the force for it to be useful; a higher value of the pre-bending angle actuation force strongly depends on the value of the gap. The peak value of the contact force is decreasing while gap is increasing.
Appl. Sci. 2019, 9, x FOR PEER REVIEW 9 of 29 Figure 6. Pre-bending angle-displacement characteristics of the V-type actuator with straight arms.
Analytical Equation (14) may be useful for the straight arms and basic design only; moreover, it underrepresents the interaction between constraint force and displacement when the displacement is close to zero or when a large force is applied [2]. From the practical point of view, actuation force should be modeled as a contact force, as in the model in [38]. Figure 7 shows the relationship between the pre-bending angle and the contact force calculated at the rigid surface, for the same model of the actuator. The gap is understood as a distance from the tip of the central shaft to the rigid surface. A small value of the pre-bending angle actuator cannot generate a high enough value of the force for it to be useful; a higher value of the pre-bending angle actuation force strongly depends on the value of the gap. The peak value of the contact force is decreasing while gap is increasing. Another important factor, which has to be considered in optimal design of the V-type actuator, is the possibility of buckling. Buckling appears if the internal force in the arms is greater than the critical buckling force. The internal compressing force in a beam is equal to the reaction force in the anchor and can be calculated from [39]: whereas critical force is given by: where min I is minor moment of inertia and β is buckling length coefficient For the verification, buckling analysis of the one actuator beam for the same FEM model was performed. Based on the same geometry, the boundary conditions for the single beam are as follows: the anchor is fixed in all degrees of freedom, while the end of the beam can move in vertical direction but cannot rotate. Critical force calculated for the geometry considered from the Equation (16) gives While solving buckling analysis by means of FEM (13), the critical load factor (CLF) was calculated. That factor is a ratio between critical and applied load for the first form of stable Another important factor, which has to be considered in optimal design of the V-type actuator, is the possibility of buckling. Buckling appears if the internal force in the arms is greater than the critical buckling force. The internal compressing force in a beam is equal to the reaction force in the anchor and can be calculated from [39]: whereas critical force is given by: where I min is minor moment of inertia and β is buckling length coefficient For the verification, buckling analysis of the one actuator beam for the same FEM model was performed. Based on the same geometry, the boundary conditions for the single beam are as follows: the anchor is fixed in all degrees of freedom, while the end of the beam can move in vertical direction but cannot rotate. Critical force calculated for the geometry considered from the Equation (16) gives F cr = 419.5 µN. While solving buckling analysis by means of FEM (13), the critical load factor (CLF) was calculated. That factor is a ratio between critical and applied load for the first form of stable model equilibrium. Assuming the same numerical data for the FEM-model simulated CLF = 1.058 indicates a good agreement between numerical and analytical results. Finally, compressive stress in the beam calculated analytically and numerically gives very similar results as well. Note, that the analytical solution incorporates only uniaxial stress and cannot provide information about local stress concentration (see Figure 2).
The calculation performed allows us to state that the numerical FEM model developed has been successfully verified. This model was used to create a parametric model for the optimization task (details in Section 5). It is apparent that the relations between geometrical parameters and the actuation force, displacement generated, buckling factor and internal stresses, are non-obvious. Moreover, for the beams with non-uniform cross-sectional area, analytical formulas cannot be applied, so the relationships between criteria are more complex.

General Information
Multiobjective optimization problems are characterized by more than one objective function. Objective functions are usually contradictory, and therefore, cannot be improved simultaneously. The objective functions should be defined in such a way that they correspond to the real demands asked of the system under consideration. In case of electrothermal microactuators it is often desired of them to have low volume, generate high force and generate low stress under working conditions, to name a few demands which can be expressed as objective functions. The determination of objective function values can be done analytically, experimentally or numerically. Analytical methods can only be applied to very simple mechanical systems (which is not the case when considering actuators of an arbitrary shape) and multiple experimental determinations of function values may be too time-consuming or even impossible, as it would require manufacturing multiple microactuators of different shapes and a complex method of measuring values related to objective functions. For those reasons, when optimizing real mechanical systems, the values of objective functions are usually determined numerically using CAE (computer aided engineering) software.
In multi-criteria problems, as opposed to single-criteria problems, where a single solution is sought, the aim is to find a set of solutions rather than a single solution. The set of solutions found should meet the design assumptions defined, as restrictions imposed on the problem and the values of the objective function should be at a level acceptable for the designer; e.g., maximal equivalent stress in the microactuator under working conditions cannot be higher than a value acceptable for a particular material used. Such an approach to the design problem is called Pareto optimization. The aim of multicriteria optimization is to find the vector of design variables describing the system, such as its geometry, material properties or boundary conditions (e.g., width of a microactuator beam, rounding radius between actuator beams and central shaft and pre-bending angle), so that the values of the functions are approaching extremes (e.g., value of vertical deflection and value of force generated). In the general case, the problem of multicriteria optimization can be presented as search of a vector of design variables.
The decision domain is, thus, n-dimensional, and the solution domain is k-dimensional. Therefore, a set of non-dominated (or Pareto-optimal) solutions is sought, also called Pareto front. Considering two vectors x and y in the decision domain, vector x strongly dominates y if: solution x weakly dominates y, if: and solution x is neutral (incomparable) to the y, if:

Visualisation of the n-Dimensional Pareto Fronts
In case of problems with two objectives, the Pareto front is a curve; in case of three objectives-a surface; and in case of higher orders of solution space-a hypersurface (Figure 8). A large number of objectives may happen to cause problems with representation of results in the graphic form, as the traditional methods may become unclear. Traditional methods of representation of result in the case of four objectives coming down to display the data as a three-dimensional plot, with three objectives shown on the axes of a coordinate system and the fourth objective displayed as a color on data points. For four objectives or more it is even harder to find an efficient way of data visualization. One way of dealing with high-dimensional results' visualization is utilizing a scatter plot matrix, where objectives are displayed pairwise, resulting in a matrix of (n 2 −n)/2 scatter plots, where n is the number of objectives. In case of six objectives, as investigated in this paper, this requires the data to be displayed on 15 scatter plots. Such an approach means that the visualization may not carried out in a clear, simple and precise manner. Utilizing self-organizing maps (SOMs) as a visualization tool helps one to overcome some of the flaws of the aforementioned methods. The decision domain is, thus, n-dimensional, and the solution domain is k-dimensional. Therefore, a set of non-dominated (or Pareto-optimal) solutions is sought, also called Pareto front. Considering two vectors x and y in the decision domain, vector x strongly dominates y if: solution x weakly dominates y, if: and solution x is neutral (incomparable) to the y, if:

Visualisation of the n-Dimensional Pareto Fronts
In case of problems with two objectives, the Pareto front is a curve; in case of three objectivesa surface; and in case of higher orders of solution space-a hypersurface (Figure 8). A large number of objectives may happen to cause problems with representation of results in the graphic form, as the traditional methods may become unclear. Traditional methods of representation of result in the case of four objectives coming down to display the data as a three-dimensional plot, with three objectives shown on the axes of a coordinate system and the fourth objective displayed as a color on data points. For four objectives or more it is even harder to find an efficient way of data visualization. One way of dealing with high-dimensional results' visualization is utilizing a scatter plot matrix, where objectives are displayed pairwise, resulting in a matrix of (n 2 -n)/2 scatter plots, where n is the number of objectives. In case of six objectives, as investigated in this paper, this requires the data to be displayed on 15 scatter plots. Such an approach means that the visualization may not carried out in a clear, simple and precise manner. Utilizing self-organizing maps (SOMs) as a visualization tool helps one to overcome some of the flaws of the aforementioned methods. SOM is an artificial neural network, which has been used, among others, in the visualization of multidimensional data sets. This method was proposed by Kohonen in [40]. In the case of multicriteria optimization problems with a large number of objectives, which result in multidimensional Pareto fronts, SOM can help to present the results in a graphical form but an alternative way. Multidimensional solution space is transformed into the form of n colorful maps, where n is the number of dimensions of solution space, and the colors of units (nodes) on the map correspond to the approximate values of the objective function. Each unit of SOM is described by a codebook vector of size n. Values of codebook vector elements are determined in an iterative learning process based on input data (Pareto front) presented to the network. The network learning process preserves the topological order of data, so solutions located close to each other on the front of Pareto remain close SOM is an artificial neural network, which has been used, among others, in the visualization of multidimensional data sets. This method was proposed by Kohonen in [40]. In the case of multi-criteria optimization problems with a large number of objectives, which result in multidimensional Pareto fronts, SOM can help to present the results in a graphical form but an alternative way. Multidimensional solution space is transformed into the form of n colorful maps, where n is the number of dimensions of solution space, and the colors of units (nodes) on the map correspond to the approximate values of the objective function. Each unit of SOM is described by a codebook vector of size n. Values of codebook vector elements are determined in an iterative learning process based on input data (Pareto front) presented to the network. The network learning process preserves the topological order of data, so solutions located close to each other on the front of Pareto remain close to each other on the SOM. Such an approach to visualization of results allows for easy determination of tradeoff relationships between objectives of the optimization problem. An effective method of visualizing a set of solutions is a key factor at the stage of information exchange between the person processing the results and the person deciding on the choice of a particular solution. In addition, SOM can be used to group similar solutions into clusters.

Formulation of the Objectives for the Optimal Design of the Microactuators
For the optimal design of the electro-thermal microactuators, the following six objectives were defined and numerically implemented: • The minimization of the volume of the structure: • The minimization of the maximal value of the equivalent stress: To calculate maximal value of the equivalent stress, the Huber-Mises-Hencky hypothesis is used [41] • The maximization of the vertical deflection of the central shaft: • The minimization of the total heat generated by the actuator: • The maximization of the buckling factor: For the calculation of the buckling factor vertical arbitrary force P v is placed to the top of the central shaft.

•
The maximization of the force generated by the actuator: The contact force is calculated with for the arbitrary value of the gap distance d g between the tip of the central shaft and the rigid surface.

Optimal Design System for Electrothermal Microactuators
In order to effectively perform multiobjective optimization of microthermal actuators, it is necessary to use appropriate optimization algorithms, define and implement optimization criteria and use efficient tools for numerical simulation. The system, which combines the immune game theory multiobjective algorithm (IMGAMO) and tools allowing for numerical simulation is proposed in the paper.
IMGAMO is the approach, which take advantages of AIS and game theory and it is used to solve the problems of multiobjective optimization [42,43].
The most important assumptions to the IMGAMO algorithm are as follows: • Each objective function (payoff function) is assigned to a single player among a group of players, playing a cooperative game, iteratively looking for a Nash equilibrium-a set of optimal solutions; • Each design variable is assigned to a single player (one player can be assigned more than one design variable); • The assignment of objective functions to players is constant, but the assignment of design variables to players changes during the course of optimization to provide a diversity of solutions; • Each player selects their strategy based on the immune algorithm to optimize its assigned objective (each player has its own population); • Each player can only modify design variables assigned to it, and the rest of design variables are fixed and determined by other players' strategies, shared between players; • The result of the algorithm (the determined Pareto frontier) is stored in the population of results.
Each objective in MOO is represented by player's payoff. Each player has its own objective (a payoff function in the Nash equilibrium). The strategy for a particular player is the optimum solution for this player's problem, remembering that other players also play their best strategies. The solution of the optimized problem consists of several parameters, each of which is assigned to one of the players. Each player optimizes only its parameters (its strategy), taking the rest of them as constant. The rest of the parameters are set according to information shared by other players. Solutions from all players should establish the solution of the problem. Then, all players use the immune algorithm to optimize their objectives.
In Figure 9, the general idea of algorithm is presented. The multiobjective optimization problem considered in the paper is solved for six criteria and seven design variables. Details of the optimization task are described in Section 5. The IMGAMO algorithm can deal with multiobjective optimization tasks for an arbitrarily chosen number of criteria and design variables. More detailed description of the IMGAMO algorithm is included in Appendix A. In Figure 9, we can see that each player is assigned one criterion and one of the decision variables (DV), which are indicated by different colors. The optimizer (clonal selection) is applied for each criterion (optimizer can change only parameters assigned to this particular player's objective). In the last step, the results are saved to a secondary population; design variables are assigned anew; and the process is repeated iteratively until the termination condition is met.
For each internal iteration in IMGAMO algorithm, consecutive evaluation of the proper criteria on the basis of encoded design variables is needed. As mentioned, particular boundary value problems have been solved by mean of FEM. MSC Mentat/Marc commercial software was adopted in order to perform those steps. The whole system of optimization consists of two main blocks: the block of the multiobjective optimization algorithm and the block of the FEM computation. Figure 10 presents the idea of the coupling between IMGAMO algorithm and block of calculation of the consecutive objectives by means of FEM in every internal iteration in IMGAMO. description of the IMGAMO algorithm is included in Appendix A. In Figure 9, we can see that each player is assigned one criterion and one of the decision variables (DV), which are indicated by different colors. The optimizer (clonal selection) is applied for each criterion (optimizer can change only parameters assigned to this particular player's objective). In the last step, the results are saved to a secondary population; design variables are assigned anew; and the process is repeated iteratively until the termination condition is met. For each internal iteration in IMGAMO algorithm, consecutive evaluation of the proper criteria on the basis of encoded design variables is needed. As mentioned, particular boundary value problems have been solved by mean of FEM. MSC Mentat/Marc commercial software was adopted in order to perform those steps. The whole system of optimization consists of two main blocks: the block of the multiobjective optimization algorithm and the block of the FEM computation. Figure 10 presents the idea of the coupling between IMGAMO algorithm and block of calculation of the consecutive objectives by means of FEM in every internal iteration in IMGAMO. Figure 10. General schema of optimal design system for electrothermal microactuators.
As mentioned in previous paragraph, B-cells of the AIS, affected by clonal selection and suppression, proliferate and generate new, possibly better solutions. Every memory B-cell contains a set of design variables. On the basis of encoded design variables, the geometry of the optimized model of the actuator is created. Next, for the built geometry, a mesh, boundary conditions, material parameters and analysis settings are established. All these steps are performed by invoking the Figure 10. General schema of optimal design system for electrothermal microactuators.
As mentioned in previous paragraph, B-cells of the AIS, affected by clonal selection and suppression, proliferate and generate new, possibly better solutions. Every memory B-cell contains a set of design variables. On the basis of encoded design variables, the geometry of the optimized model of the actuator is created. Next, for the built geometry, a mesh, boundary conditions, material parameters and analysis settings are established. All these steps are performed by invoking the Mentat preprocessor.
Calculation of the particular k-th objective is performed on the basis of output file generated by MSC Marc. It should be underlined that solution times for particular analyses differ between each other. As in the proposed method of optimization, each criterion (represented by player's payoff) is optimized independently; the optimization procedure can be much more efficient compared to the other optimization methods, for which in every iteration all the values of objective functions have to be determined for each solution.

Optimization Model
The model of the structure which contains one pair of the arms was considered for the shape optimization. Symmetry along the vertical axis in the center shaft was assumed, which ensured axial movement of the central shaft. The beam was modeled by means of the NURBS curve, which consists of four control points. The knot vector of the NURBS curve is uniformly distributed and has cubic interpolation; the constant value of the control points' weight being equal to 1.0 is assumed. The shape of the beam is modeled by manipulation of the NURBS control polygon, with symmetry imposed along the centerline of the beam. The first four design variables are distances between a particular control point and the symmetry line. The fifth, sixth and seventh design variables are the pre-bending angle of the beam, radius R1 and radius R2, respectively, so the total number of design variables is seven. Figure 11 presents geometry of the model and the way of parametrization. Calculation of the particular k-th objective is performed on the basis of output file generated by MSC Marc. It should be underlined that solution times for particular analyses differ between each other. As in the proposed method of optimization, each criterion (represented by player's payoff) is optimized independently; the optimization procedure can be much more efficient compared to the other optimization methods, for which in every iteration all the values of objective functions have to be determined for each solution.

Optimization Model
The model of the structure which contains one pair of the arms was considered for the shape optimization. Symmetry along the vertical axis in the center shaft was assumed, which ensured axial movement of the central shaft. The beam was modeled by means of the NURBS curve, which consists of four control points. The knot vector of the NURBS curve is uniformly distributed and has cubic interpolation; the constant value of the control points' weight being equal to 1.0 is assumed. The shape of the beam is modeled by manipulation of the NURBS control polygon, with symmetry imposed along the centerline of the beam. The first four design variables are distances between a particular control point and the symmetry line. The fifth, sixth and seventh design variables are the pre-bending angle of the beam, radius R1 and radius R2, respectively, so the total number of design variables is seven. Figure 11 presents geometry of the model and the way of parametrization. The actuator was modeled as plane stress structure fabricated from polycrystalline silicon (Young's modulus E = 158 × 10 3 MPa, Poisson's ratio ν = 0.23, thermal conductivity k = 140 × 10 8 pW/µmK, resistivity R = 3.3 × 10 −11 TΩµm and linear expansion coefficient α t = 3 × 10 −6 1/K). The length of the beam was 220 µm, whereas the thickness of the in the z-direction was 10 µm (µMKS system unit is used). The model was supplemented with electrical, thermal and mechanical boundary conditions. Anchors were fixed in all degrees of freedom and the temperature 300 K was assumed. The electric potential difference between anchors was equal to 5 V. For the buckling analysis force P v = 500 µN, whereas for the calculation of the contact force, the gap d g was equal to 2 µm. Box constraints imposed on the design variables were as follows: for the design variables D1, D2, D3, D4, R1 and R2 within the scope < 1 ÷ 20 > [µm], whereas the fifth design variable ϕ was within the scope, < 0.001 ÷ 12 > [ • ]. The multiobjective optimization task concerned determining seven design variables which minimized six defined functionals (Section 3.3) and was performed by means of IMGAMO algorithm. A few runs of IMGAMO algorithm were performed. Due to the simplicity, only two are considered and described in the paper. The following values of IMGAMO parameters were assumed:

Results and Discussion
For the first optimization task (variant 1), a six-dimensional Pareto front consisting of 2471 nondominated solutions was found, whereas for the variant 2, there were 3143 solutions. Figures 12  and 13 present scatter plot matrices generated on the basis of set of Pareto optimal solutions, including the distribution of objectives' values in the form of histograms. To transform maximization problems into minimization problems, some of them were multiplied by −1 so that only minimization tasks were solved.      The size of the set of Pareto optimal solutions is different for both variants; however, shapes of the Pareto fronts in each projection plane are very similar. That shows that it is highly probable, that the solutions obtained are unambiguous and are located close to the optimal front. The differences can only be seen in the filling of the space by the nondominated solutions. It has to be underlined, that the computational effort for the first variant is significantly smaller. A detailed comparison can be found in the Table 1. Table contains a   In order to obtain the scatter plot matrix present, by means of popular bi-objective optimizers, 15 independent runs of optimization tasks were needed. Moreover, for that optimization problem, a multidimensional Pareto front could not be obtained by solving several separate bi-objective tasks, because in this case, structures with no contact force are ignored, and they would not be included in the set of nondominated solutions. Due to this, some comparisons are possible only for the projection plane, which incorporates objective f6. In attempt to compare the results of IMGAMO algorithm with those of a popular and well-known optimization algorithm, a bi-objective NSGAII was chosen as a reference algorithm. A pair of f2 and f6 objectives was selected for the comparison. Figure 14 shows Pareto front obtained by NSGAII, and the projection plane for criteria f2 and f6, obtained by method proposed herein (such projection plane is indicated in Figure 13). The NSGAII algorithm was run for the population size of 50 chromosomes and 100 iterations. (File containing data with Pareto optimal solutions in each generation are included in supplementary material). As can be seen, the location of the Pareto front found by both algorithms is comparable; however, IMGAMO found more a widespread one. This comparison also shows the high probability of the correct location of the Pareto front for this problem, found using the method proposed herein. because in this case, structures with no contact force are ignored, and they would not be included in the set of nondominated solutions. Due to this, some comparisons are possible only for the projection plane, which incorporates objective f6. In attempt to compare the results of IMGAMO algorithm with those of a popular and well-known optimization algorithm, a bi-objective NSGAII was chosen as a reference algorithm. A pair of f2 and f6 objectives was selected for the comparison. Figure 14 shows Pareto front obtained by NSGAII, and the projection plane for criteria f2 and f6, obtained by method proposed herein (such projection plane is indicated in Figure 13). The NSGAII algorithm was run for the population size of 50 chromosomes and 100 iterations. As can be seen, the location of the Pareto front found by both algorithms is comparable; however, IMGAMO found more a widespread one. This comparison also shows the high probability of the correct location of the Pareto front for this problem, found using the method proposed herein. Extreme solutions for each objective functional are collected in Table 2, whereas Figure 15 presents ghd optimal shape of the actuator for that extremes. Additionally, for comparison, actuators with straight arms (traditional geometry) for different widths of the beam equal to 2.5 μm , 5.0 μm Figure 14. Comparison between Pareto front obtained by NSGAII projection plane for criteria f2 and f6 obtained by means of IMGAMO algorithm.
Extreme solutions for each objective functional are collected in Table 2, whereas Figure 15 presents ghd optimal shape of the actuator for that extremes. Additionally, for comparison, actuators with straight arms (traditional geometry) for different widths of the beam equal to 2.5 µm, 5.0 µm and 10.0 µm are included. For such geometry, the pre-bending angle was assumed to ϕ = 1.3 • and a sharp connection between arm and the anchor was assumed. It is natural, that a significantly higher maximal value of the equivalent von Mises stress was observed. Table 3 contains the value of the design variables obtained for the extreme solutions. Representation of the set of Pareto optimal solutions by means of scatter plot matrix, provides good insight into the nature of the problem. Shapes of the consecutive tradeoffs show dependences between particular objectives. The rapid change at some region of the Pareto front can be very useful information for the designer (small changes in one criterion may cause great changes in the other). However, a higher number of criteria, analyzing too many charts simultaneously, can be a serious obstacle in the effective selection of the interesting solutions. As previously mentioned, utilizing the SOMs as a visualization tool may help to overcome that difficulty and could be very useful for the designer.  Representation of the set of Pareto optimal solutions by means of scatter plot matrix, provides good insight into the nature of the problem. Shapes of the consecutive tradeoffs show dependences between particular objectives. The rapid change at some region of the Pareto front can be very useful information for the designer (small changes in one criterion may cause great changes in the other). However, a higher number of criteria, analyzing too many charts simultaneously, can be a serious obstacle in the effective selection of the interesting solutions. As previously mentioned, utilizing the SOMs as a visualization tool may help to overcome that difficulty and could be very useful for the designer.
On the basis of set of Pareto optimal solutions, colorful grid maps by means of SOM were generated. Figure 16 presents the maps generated for the 2nd variant of optimization. For the 1st variant, the maps generated are comparable, so they are not included in the paper. Each grid in the  On the basis of set of Pareto optimal solutions, colorful grid maps by means of SOM were generated. Figure 16 presents the maps generated for the 2nd variant of optimization. For the 1st variant, the maps generated are comparable, so they are not included in the paper. Each grid in the figure corresponds to a particular optimization criterion, whereas nodes correspond to particular solutions, with the value of each objective depicted as the unit's color. The general look and the brief analysis of the maps allow us to formulate conclusions:

•
Objectives f1 and f4 have quite a similar nature (actuator with a smaller volume generates much smaller amounts of heat); • Objectives f3 and f5 have contradictory natures, except in the region in the upper right corner of the grid (geometry of the arms resistant to buckle can generate low displacement and vice versa), • Objectives f4 and f5 have contradictory natures as well.
Deeper insight into the particular regions of the SOM makes it possible to identify more detailed dependencies between the criteria. Now consider the process of selecting interesting compromise-solutions on the basis the SOM generated. Let us start from the analysis of the criterion f2, by rejecting the solutions with higher values of the maximal equivalent stress (step 1 in Figure 17). The average tensile strength of polysilicon is about 1.45 GPa, so the assumed threshold level was set to 1000 MPa. Then, we got rid of the solution which could buckle. The threshold level was set to the buckling factor about value 2.0. Transmitting such limited region to the criterion f3 (step 2), the two regions with higher value of generated displacement by the actuator can be identified. For the region located in the center (heart shaped curve), smaller values of criteria f1 and f4 were observed, compared to the second region (oval shaped). If we accept the level of the value of actuation force (objective f6), the solution from that region may be chosen (step 3). Figure 18 shows the geometry of the actuator representing that region, whereas values of the objectives and design variables for that solution are included in Tables 2 and 3. generated displacement by the actuator can be identified. For the region located in the center (heart shaped curve), smaller values of criteria f1 and f4 were observed, compared to the second region (oval shaped). If we accept the level of the value of actuation force (objective f6), the solution from that region may be chosen (step 3). Figure 18 shows the geometry of the actuator representing that region, whereas values of the objectives and design variables for that solution are included in Tables 2 and 3.   Note that this compromise-solution is similar to the shape of the actuator proposed by Kwan A. et al., in the work dedicated to improving designs for an electrothermal in-plane microactuator [7]. Note that this compromise-solution is similar to the shape of the actuator proposed by Kwan A. et al., in the work dedicated to improving designs for an electrothermal in-plane microactuator [7]. The actuators with nonuniform beams, for which width increases close to the central shaft, allow for the generation higher of displacement and actuation force compared to the actuator with uniform beams. Other examples of selecting compromise-solutions with respect to the different assumption and restrictions, starting from another arbitrarily chosen criterion, can be found as well.  Note that this compromise-solution is similar to the shape of the actuator proposed by Kwan A. et al., in the work dedicated to improving designs for an electrothermal in-plane microactuator [7]. The actuators with nonuniform beams, for which width increases close to the central shaft, allow for the generation higher of displacement and actuation force compared to the actuator with uniform beams. Other examples of selecting compromise-solutions with respect to the different assumption and restrictions, starting from another arbitrarily chosen criterion, can be found as well.
It should emphasized, that solutions for which the maximal equivalent stress exceeds strength of the material, and solutions for which buckling may occur (factor less than 1.0), are infeasible. That may lead to the conclusion that such criteria should be implemented in optimization as constraints. It is of course possible and can be easily implemented in the method of optimization proposed, and would cause a reduction of the total number of objectives. On the other hand, as previously mentioned, making no assumptions on weights of objectives before optimization, can result in obtaining useful information for the designer concerning trade-off relationships between criteria. Moreover, there are no criteria related to the dynamic characteristic of the actuator considered in this paper. The method of multicriteria optimization proposed allows for the formulation and implementation of more criteria as well.

Final Conclusions
A novel method for the multiobjective optimization of an electrothermal microactuator has been presented. A novelty of the proposed method is particularly manifested in effective shape optimization related to the large number of criteria. Six different criteria were formulated and It should emphasized, that solutions for which the maximal equivalent stress exceeds strength of the material, and solutions for which buckling may occur (factor less than 1.0), are infeasible. That may lead to the conclusion that such criteria should be implemented in optimization as constraints. It is of course possible and can be easily implemented in the method of optimization proposed, and would cause a reduction of the total number of objectives. On the other hand, as previously mentioned, making no assumptions on weights of objectives before optimization, can result in obtaining useful information for the designer concerning trade-off relationships between criteria. Moreover, there are no criteria related to the dynamic characteristic of the actuator considered in this paper. The method of multicriteria optimization proposed allows for the formulation and implementation of more criteria as well.

Final Conclusions
A novel method for the multiobjective optimization of an electrothermal microactuator has been presented. A novelty of the proposed method is particularly manifested in effective shape optimization related to the large number of criteria. Six different criteria were formulated and numerically implemented. Multiobjective shape optimization for such criteria was successfully solved. New, additional objectives can be formulated and implemented for the method of optimization proposed. One of the main advantages of the method proposed is the possibility of obtaining a multidimentional set of Pareto-optimal solutions by performing a single optimization task. The IMGAMO algorithm significantly outperforms other popular multiobjective optimization algorithms (NSGAII and SPEA2), when a large number of criteria is considered. Several numerical tests on the mathematical benchmark test functions have been performed. A detailed comparison is included in the paper.
From a practical point of view, solving the problem considered by means of scalarization method, ε-constrained method or performing several separate bi-objective problems, can be inadequate and problematic. The superiority of method proposed here in comparison to the other approach for the problem considered has been presented and discussed. The method proposed for the parametrization of the actuator's beam allows a reduction in the total number of design parameters and allows for generating shapes with significantly lower stress concentrations. The multidimensional Pareto fronts obtained may consist of thousands of solutions. Selecting the proper solution by browsing such a big data set is very difficult and inefficient. A representation of the multidimentional Pareto fronts by means of scatter plot matrices and aided by SOM has been presented. The method of selecting interesting compromise-solutions has been proposed and described. It can be very helpful as an alternative tool in the optimal design of microthermal actuators. The method proposed for the optimization of a large number of, can be used for solving other multiobjective and multidisciplinary problems related to the different engineering disciplines.
It can be a very time-consuming task, especially for a more complicated geometry, when more design variables or more objectives are considered. Moreover, the solution of the coupled problems is more time-consuming when compared to a single-field problem. Parallel computation should be the next step in order to reduce the time of the optimization.

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

Appendix A
This section contains additional information related to the details of the IMGAMO algorithm. The algorithm works for an arbitrary (determined by user) number of design variables and number of objectives. To run the algorithm, a few parameters have to be set. The first of them is population size (popSize)-the size of the population for each player is stored. Another parameter is the number of clones (clonNumb)-number of clones which will be created in each iteration. Iterations store data about number of iterations in algorithm, and dSup-is suppression distance. It is the minimal distance between individuals in objective space. Such parameters should be selected, taking into account the difficulty of the optimization problem. Problems with higher numbers of design variables and higher numbers of objectives are understood to be more difficult. Very often, practical engineering problems belong to the class of NP-hard optimization problems, so the application of the metaheuristic approaches can be the only solution in this situation. Determination of the values of the parameters allows for control of the exploitation near the probable optimum and the exploration of the search space. Moreover, the user has an influence of the size of the set of optimal Pareto solutions. It should be underlined that, due to the probabilistic character of the IMGAMO algorithm (like for all metaheuristic approaches), different runs with the same value of the parameters do not guarantee obtaining the same results. Figure A1 presents the schema of the IMGAMO algorithm. IMGAMO algorithm sequence of operations, which is graphically presented in Figure 10 can be expressed as follows: k is the number of objectives in the problem; 2.
k players are created and populations initialized with size population size are created for each of them; Figure A1. General schema of IMGAMO algorithm.
IMGAMO algorithm sequence of operations, which is graphically presented in Figure 10 can be expressed as follows: 1.
k is the number of objectives in the problem; 2.
k players are created and populations initialized with size population size are created for each of them; 3.
Empty population result is created; 4.
Each design variable of the problem is randomly assigned to one of the players; 5.
Each player runs the clonal selection procedure: a.
Calculate the fitness (assigned t player) of the solution; b.
Sort the population taking into account the fitness; c.
The best solution creates clon_number clones; d.
Each next solution creates one clone less; e.
Each clone is mutated (only parameters optimized by this player are mutated)-hipersomatic mutation f.
If the fltness of the mutated clone is better than the parent replace parent with clone.

9.
If exchange iter = clon_iter, a. Take the best solutions from each player; b.
Store optimized parameters from each player to other populations as constant; c.
Each parameter of the problem is randomly assigned to one of the players.

Appendix B
This section presents the effectiveness of the IMGAMO algorithm. Several tests were performed to analyze and check the effectiveness of IMGAMO algorithm with respect to one of the most popular metaheuristic optimization algorithms (NSGA II and SPEA2). For the comparison, the set of DTLZ test functions were chosen [44]. Formulations of the DTLZ problems are presented in Table A1. SPEA2 and NSGAII approaches with parameters in Table A2. Each algorithm was run 10 times and the average values of metrices were calculated. In order to evaluate algorithm, two metrices were chosen: • Generational Distance (GD): • Spacing (S): where Apart from those two metrics, the comparison was also performed for the size of the set of Pareto-optimal solutions (Frontsize).
For the DTLZ1 problem (Table A3), one can find out that the IMGAMO algorithm significantly outperformed its competitors. There are a very low values for GD and S metrices. For the DTLZ2 problem (Table A4), the IMGAMO algorithm also outperformed the others. For six-criteria DTLZ2 experiments with lower numbers of functions, evaluations were performed. One can find that even for 1000 evaluations, results are satisfactory and better than NSGAII and SPEA2 with 100,000 evaluations (Table A5). This four-criteria problem was also formulated with 50 decision variables. Again, for that experiment, IMGAMO gave better results. For the next problem-DTLZ3-even for a three criteria formulation, IMGAMO significantly outperformed NSGAII and SPEA2 (Table A6). For the experiment with higher number of criteria, the results of IMGAMO were also satisfactory. The next problem considered-DTLZ4, for all formulations, lead to similar results. IMGAMO gave the best results for all experiments (Table A7). For the last problem analyzed-DTLZ5-the experiments were performed for four, five, six and seven criteria in comparison to SPEA2 and NSGA-II, and additionally, for eight, nine and 10-criteria formulations (Table A8). All experiments proved that the IMGAMO algorithm deals with such problems, finding fronts close to the real ones. To summarize, the biggest strength of this approach is in dealing with problems with higher numbers of criteria. Algorithm IMGAMO has successfully found the fronts close to the real ones for problems with 4-10 criteria, dealing with several difficulties:
For some cases, IMGAMO found a less numerous set of the Pareto optimal solutions (FrontSize), but this factor has significant importance if solutions are close to the true Pareto front (GD) and are well distributed (S). For the cases where IMGAMO found smaller Frontsizes than NSGA II and SPEA2, values of metrices GD and S were significantly worse. That means that much more numerous FrontSize could be a local front or solutions are not well distributed.
The algorithm showed its scalability in decision variables' space. In problems DTLZ2 and DTLZ4, the number of variables was increased to 50. The IMGAMO algorithm, without any problems, found the front sought.