## 5 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. Figure 11: Geometry and parametrization of the model. The actuator was modeled as plane stress structure fabricated from polycrystalline silicon (Young's modulus \(\mathrm{E}=158\,\times\,10^{3}\) MPa, [PERSON]'s ratio \(\
u=0.23\), thermal conductivity \(\mathrm{k}=140\times 10^{8}\)\(\mathrm{p}\mathrm{W}/\mathrm{\SIUnitSymbolMicro}\mathrm{m}\mathrm{K}\), resistivity \(\mathrm{R}=3.3\times 10^{-11}\mathrm{T}\mathrm{\SIUnitSymbolMicro}\mathrm{m}\) and linear expansion coefficient \(\alpha_{t}=3\times 10^{-6}\,1/\mathrm{K}\)). The length of the beam was 220, whereas the thickness of the in the z-direction was 10 (, whereas 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\), whereas for the calculation of the contact force, the gap \(d_{g}\) was equal to 2. 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\div 20>[\mathrm{\SIUnitSymbolMicro}\mathrm{m}]\), whereas the fifth design variable was within the scope, \(<0.001\div 12>[\mathrm{\SIUnitSymbolMicro}\mathrm{\SIUnitSymbolMicro}\mathrm{ \SIUnitSymbolMicro}  
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 number of fitness function evaluations for each objective (columns 2-7), the total number of fineness function evaluations (column 8) and the size of the set of Pareto optimal solutions (column 9). Files containing data (values od design variables and values of objective functions) from the set of Pareto-optimal solutions are included in supplementary material. 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.