Optimization of Generatively Encoded Multi-Material Lattice Structures for Desired Deformation Behavior

: Natural systems achieve favorable mechanical properties through coupling signiﬁcantly different elastic moduli within a single tissue. However, when it comes to man-made materials and structures, there are a lack of methods which enable production of artifacts inspired by these phenomena. In this study, a method for design automation based on alternate deposition of soft and stiff struts within a multi-material 3D lattice structure with desired deformation behavior is proposed. These structures, once external forces are applied, conform to the geometry given in advance. For that purpose, a population-based algorithm was proposed and integrated with a multi-material physics simulator. To reduce the amount of data processed during optimization, a generative encoding method based on discrete cosine transform (DCT) was proposed. This enabled a compressed topological description and promoted symmetry in material distribution. The simulation results showed different three-dimensional lattice structures designed with proposed algorithm to meet a set of desired deformation behaviors. The relation between residual deformation error, targeted deformation geometry, and material distribution is discussed.


Introduction
Processes of manufacturing and design are tightly coupled. The approach designers take when shaping a product is biased by both CAD (Computer Aided Design) platform used, but also by manufacturing process necessary for product realization. Usually, so called feature-based modeling, in which either standard 2D or 3D primitives are used to create 3D volumes out of homogenous material by standard CAD tools, are applied. These standard CAD environments were not created with a freeform design paradigm in mind. This is a limiting factor for their accommodation to complex, non-homogenous, multimaterial 3D design and analysis [1,2]. The standard design approach has its justification in manufacturing costs, only if the technologies based on additive manufacturing, that naturally create complex freeform shapes of arbitrarily internal material distribution, are omitted from consideration. These latter technologies can deliver products of the highest complexity, where the complexity of the manufactured part comes at no added cost. This opens a vast space for engineers to explore and utilize when designing completely new products, with new and previously unattainable properties. Lattice structures are especially suited for manufacturing with additive technologies. Their complex internal morphology is materialized through layer-by-layer material deposition, which is hard to replicate with conventional manufacturing processes. Lattices have several advantages over homogenous structures, due to their favorable mass to strength ratios. These complex structures have received more attention recently in the scientific community because they can be used to fill out the internal space of a volume defined by arbitrarily complex freeform surfaces [3].
A considerable amount of work has been done in the field of topology and shape optimization, especially for mono-phase material structures. Both lattice-based and continuous volumes were considered, with the most common optimization criteria of compliance ing the analyzed soft and stiff materials, we show in the following sections how certain combinations of materials limit the ability of the structure to comply to a desired deflection line. The structure shown in Figure 1a is a homogenous lattice structure, with illustrated loads and restraints. It deforms in an expected way, referred to as the normal response. On the right side is a structure optimized by deposition of soft and stiff lattice struts with the goal of deforming in a straight line. Theoretically, an indefinite number of phases exist through continuous mixing of resins, but only a selected combination can be printed in a single structure, due to 3D printers' manufacturers limitations, and physical compatibility of distinct phases. Regarding the analyzed soft and stiff materials, we show in the following sections how certain combinations of materials limit the ability of the structure to comply to a desired deflection line. The structure shown in Figure 1a is a homogenous lattice structure, with illustrated loads and restraints. It deforms in an expected way, referred to as the normal response. On the right side is a structure optimized by deposition of soft and stiff lattice struts with the goal of deforming in a straight line. The deflection lines considered in the paper range from simple ones such as straight line, to complexly composed ones including combinations of linear members, negative curvatures, step deflection, and fourth order polynomials, where locations of inflections and slopes are manipulated. These unintuitive responses have the potential to be used in different applications, i.e., ductility enhancements, fracture toughness, negative to zero thermal expansion, biomimetic materials, soft robotics, and functional materials [19][20][21][22]. Comparable concepts are also present in biological tissues, bones, plants, etc., where nature has produced a strategy to achieve unusual mechanical properties through coupling variable elastic moduli from a few GPa to below KPa [23].
The problem of finding material distribution to meet desired profile is formulated in optimization domain. A robust optimization algorithm is needed to reliably explore this domain. An evolutionary algorithm is developed and tested for the purpose in this study. Since the problem considered here has a vast search landscape, which is related to the number of struts in lattice and number of phases of materials used, a particular consideration in the development of the algorithm was dedicated to the genotype-phenotype mapping scheme since it is a time-consuming procedure.
It is necessary to find such an encoding that will allow the representation of the whole structure using a minimal number of parameters. The scheme should be simple, but allow for the representation of a desired level of details in the structure. It should also naturally promote evolvability of the population. For that reason, a generative encoding procedure was proposed, based on inverse Discrete Cosine Transform (DCT).
This procedure compresses the data necessary to describe complete structures at the arbitrarily fine resolution, and has some other favorable traits, such as natural symmetry promotion, which is important since it allows reduced calculation of error profiles, as explained in detail in Section 2.2.
The main contributions in this paper can be summarized as follows: (1) Proposing a method based on the Discrete Fourier Transform (DFT) for compressing the topological data of multi-material lattice structures, and reliable transformation from genotypic to phenotypic space and vice versa. This allows reduction of the convergence time in the evolutionary optimization domain. The method described in this study can be considered as a generalization of methods proposed to describe continuous topology in the frequency The deflection lines considered in the paper range from simple ones such as straight line, to complexly composed ones including combinations of linear members, negative curvatures, step deflection, and fourth order polynomials, where locations of inflections and slopes are manipulated. These unintuitive responses have the potential to be used in different applications, i.e., ductility enhancements, fracture toughness, negative to zero thermal expansion, biomimetic materials, soft robotics, and functional materials [19][20][21][22]. Comparable concepts are also present in biological tissues, bones, plants, etc., where nature has produced a strategy to achieve unusual mechanical properties through coupling variable elastic moduli from a few GPa to below KPa [23].
The problem of finding material distribution to meet desired profile is formulated in optimization domain. A robust optimization algorithm is needed to reliably explore this domain. An evolutionary algorithm is developed and tested for the purpose in this study. Since the problem considered here has a vast search landscape, which is related to the number of struts in lattice and number of phases of materials used, a particular consideration in the development of the algorithm was dedicated to the genotype-phenotype mapping scheme since it is a time-consuming procedure.
It is necessary to find such an encoding that will allow the representation of the whole structure using a minimal number of parameters. The scheme should be simple, but allow for the representation of a desired level of details in the structure. It should also naturally promote evolvability of the population. For that reason, a generative encoding procedure was proposed, based on inverse Discrete Cosine Transform (DCT).
This procedure compresses the data necessary to describe complete structures at the arbitrarily fine resolution, and has some other favorable traits, such as natural symmetry promotion, which is important since it allows reduced calculation of error profiles, as explained in detail in Section 2.2.
The main contributions in this paper can be summarized as follows: (1) Proposing a method based on the Discrete Fourier Transform (DFT) for compressing the topological data of multi-material lattice structures, and reliable transformation from genotypic to phenotypic space and vice versa. This allows reduction of the convergence time in the evolutionary optimization domain. The method described in this study can be considered as a generalization of methods proposed to describe continuous topology in the frequency domain [24,25]. (2) Multi-material lattice domain space exploration by evolutionary algorithms allows optimization of the complete structures' deformation behavior, as opposed to deformation performed in a single point.
The whole lattice structure with all attributes, including topology, elasticity constants, and shape, is uniquely mathematically described using connectivity and nodal matrices, is scalable in terms of including additional materials, and is suitable for manufacturing using additive manufacturing processes.
The remainder of the paper is organized as follows: the methods and algorithm outline are presented in Sections 2 and 3. In Section 4, a number of 3D structures optimized to conform to deflection lines of different complexities are presented. Limitation of the structures based on material properties and complexity of deflection lines are discussed. In Section 5, discussion is given, and results are compared to similar existing approaches. Section 6 gives conclusions and insights in future research.

Evolutionary Algorithms
Evolutionary algorithms belong to the family of population based heuristic optimization algorithms. These algorithms are inspired by some of the processes found in biological evolution, i.e., selection, mutation, and recombination. These mechanisms are used to evolve, or to artificially breed a solution for a typically hard optimization problem. The basic principle here is that individuals that form a population of potential solutions compete to survive. An individuals' fitness, which is proportionate to its survival rate, is measured with the fitness function. Fitter individuals have a higher chance to survive and to pass their genetic material through crossover and mutation to the following generations.
Evolutionary algorithms, and other population-based algorithms, have been successfully applied to different problem domains. Nondeterministic polynomial time (NP) hard problems, large search spaces, constrained optimization problems with deceptive domains, among others, have all been successfully explored and optimized with evolutionary algorithms [26][27][28]. The disadvantages of evolutionary algorithms (EA)s compared to other optimization techniques lie in their stochastic nature. There is no formal guarantee of finding global optima. The quality of the solution depends on many factors, whose interconnections are highly nonlinear and stochastic in nature. Thus, finding a satisfactory combination of these parameters, i.e., population size, selection and surviving mechanism, probabilities of mutation and recombination, fitness scaling, genotype to phenotype mapping, etc., is very problem dependent and requires tedious experimental fine tuning to yield satisfactory results. Additionally, similarly to the biological process of evolution, the evolutionary algorithms might become slow for larger populations of individuals holding large information size in their genotype.
In this study, we propose a procedure based on generative encoding concept to compress the data contained in each individual, using a frequency matrix, and an inverse discrete cosine transform to make the mapping from genotype-this is a computational domain where crossover and mutation take place, to phenotype-this is the physical space where selection takes place. In our example, genotypic space consists of 3D frequency matrices where frequencies dictate material properties of struts, while phenotypic space consists of actual three-dimensional multi-material lattice structures subjected to loads and restraints.

Generative Encoding
It is worth describing encoding procedure in more detail. In generative encodingbased approaches, the information about the individual in the population are stored indirectly applying a set of rules. These rules are used to perform the transition of the individual from the phenotypic to genotypic space and vice versa. This can be considered a way of compressing the amount of data that describes the individual. The drawback of this approach is that there exists a computational overhead in applying these rules. Some examples of successfully applying generative encodings are generative grammar approaches, which were used to reduce the search space in evolving conceptual tensegrity structure deigns [29]. Compositional pattern-producing networks (CPPN) were used to optimize shape and material distribution of free-form objects [30]. Gaussian Mixtures Symmetry 2021, 13, 293 5 of 13 (GMX) and its derivatives, based on level-set approaches, have been used to describe density fields and successfully applied to describe continuous soft-body robot systems [31].
In this study, we use a DCT approach to represent material distribution of the lattice structure. It is applied here for its simplicity and proven evolvability, and the ability to easily describe increasingly complex material distributions by including higher order harmonics. The distribution of material along a specific axis can be interpreted as a onedimensional signal and treated as Fourier series. The usefulness of the Fourier analysis is that we can break up any arbitrary periodic function into a set of simple terms that can be solved individually and then recombined to reconstruct the original signal to a high degree of accuracy.
In the approach adopted in this paper, a phenotype of an individual is a 3D matrix composed of frequency amplitudes ranging from an arbitrarily chosen interval. Each amplitude is related to an adjacent harmonic, and the number of harmonics for each axis included in the phenotype limits the complexity of the design that can be reconstructed from the frequency matrix for that specific axis.
If there are N harmonics dedicated to the specific axis of the lattice design, there are potentially N regions along this axis, each describing a region of different mechanical properties of the structure. If we sum the numbers of all harmonics for all axes, we can define a metric called fineness of the structure: where j is for all the harmonics associated with a particular axis, and i is for the number of axes considered: 1 is for 1D case, 2 is for 2D, and 3 is for a 3D case. The Discrete Cosine Transformation sums all weighted sinusoids and applies a threshold at a given value to define distinct regions of materials. This procedure is applied for each axis of the structure allowing the reconstruction of volumetric density field for the complete lattice structure. The structure illustrated in Figure 2 is represented by summation of three weighted sinusoids to form a three-dimensional weighted space. The three sine waves have the following form: Symmetry 2021, 13, x FOR PEER REVIEW 6 of 13 xi denote the midpoint of the ith segment and let fi denote f (xi). To avoid aliasing, the values of k are limited to be k = 0, 1, 2, …, N-1. The integral (4) is approximated by the midpoint quadrature rule, yielding: which is, in some settings, known as Discrete Cosine Transform of Type II (DCT-II). This is the most used form and referred to as the DCT [32]. It is important to stress here that describing the structure using the approach presented in this study not only compresses the data required for representation of the structure, but also promotes symmetry of material deposited through the structure. This is exploited for the case of symmetrical loads and restraints analyzed in this study. In general, to meet the targeted defection profile, at least four contours, each associated to one edge along the longest dimension, in our case, x1 of the structure, should be evaluated. This also holds for the case of symmetrical loads and restraints, but without promoting the symmetry through appropriate encoding. In the approach presented in this paper, all the deflection error calculations were performed only for one edge, which proved enough, as will be explained in the following sections.
Before rendering a multi-material lattice structure, it is necessary to check the position of the center of mass (cm) of each strut. The position of cm for the strut i in the 3D This space is thresholded by the usage of inverse discrete cosine transformation. In the material distribution given above, the lattices belonging to regions above the threshold are defined as stiff and vice versa. This approach easily represents 3D topologies with arbitrarily detailed resolution of geometry in all three principal axes. The one-dimensional Fourier cosine series on the interval [0, L] is given by: with coefficients: This representation is useful for even functions. Now, the discrete version of the Fourier cosine series can be derived by dividing the interval [0, L] into N equal segments. Let x i denote the midpoint of the ith segment and let f i denote f (x i ). To avoid aliasing, the values of k are limited to be k = 0, 1, 2, . . . , N-1. The integral (4) is approximated by the midpoint quadrature rule, yielding: which is, in some settings, known as Discrete Cosine Transform of Type II (DCT-II). This is the most used form and referred to as the DCT [32]. It is important to stress here that describing the structure using the approach presented in this study not only compresses the data required for representation of the structure, but also promotes symmetry of material deposited through the structure. This is exploited for the case of symmetrical loads and restraints analyzed in this study.
In general, to meet the targeted defection profile, at least four contours, each associated to one edge along the longest dimension, in our case, x 1 of the structure, should be evaluated. This also holds for the case of symmetrical loads and restraints, but without promoting the symmetry through appropriate encoding. In the approach presented in this paper, all the deflection error calculations were performed only for one edge, which proved enough, as will be explained in the following sections.
Before rendering a multi-material lattice structure, it is necessary to check the position of the center of mass (c m ) of each strut. The position of c m for the strut i in the 3D density matrix defines the property of the material of this strut, proportional to the material density field. This procedure must be repeated for each strut in the lattice structure. Despite that, the computational overhead in doing this is not excessively intensive, because the topology of the lattice stays fixed over time, since the nodes remain fixed and there is no variation of the length of the struts. For that reason, no iteration over the structure is necessary, a lookup table is designed in advance holding the information of c m and compared to mass distribution matrix.

Fitness Evaluation
Fitness function presents a central part of any evolutionary algorithm. This function is, in evolutionary context, an environment to which competing individuals must adapt. The ones that can make this adaptation better are fitter and have an increased chance to survive compared to the individuals that are less fit. In the example presented here, individuals must evolve a material distribution which enables them to achieve desired structural deformation.
To measure how a multi-material lattice structure deforms, the following is needed: a method to create the outer shape of the structure, a method to fill this shape with lattice structure, and a method to apply restraints, loads, and evaluate structural response.
As it concerns the first problem, outer geometry defined as a STEP model method is applied in this paper. STEP-file can support arbitrarily complex shapes, is free, easy to read, and recognized through the ISO-10303-21 standard. The outer geometry then is used to bound the lattice. The lattice itself is constructed using a standard tetrahedral mesh generator. In this study Gmsh generator [33] is used.
The output from the mesh generator consists of two matrices: connectivity matrix, later referred to as Conn.txt, and a nodal coordinate matrix, later referred to as Coord.txt. These two matrices are fed to the next step. In the next step, material and volume is associated with each strut. This enables the calculation of the response of the structure once loads are applied. The following step is to employ a multi-material physics simulator to evaluate structures' response based on the distribution of the material, loads and restraints. Standard nonlinear finite element solvers are of limited applicability in situations including large deformations, instabilities, and significant differences between elastic moduli of struts [34]. Another problem is their computational time, which makes them impractical for implementation in population-based optimization algorithms. For that purpose, Vox-CAD [35][36][37][38][39][40], a multi-material physics simulator is used. Voxelyze, the underlying engine running behind VoxCAD, is a computationally efficient simulation approach based on nonlinear relaxation of structural elements. Voxelyze simulates elastic voxels based on an internal lattice of discrete points, interconnected by spring like beam elements with translational and rotational stiffness. Once it is possible to calculate deformation of the structure subjected to loads, the optimization algorithm distributes the soft and stiff lattice struts throughout the structure, trying to minimize the deviation between the actual deflection profile, and the desired deflection profile. The outline of the optimization processes are illustrated in Figure 3. a method to create the outer shape of the structure, a method to fill this shape with lattice structure, and a method to apply restraints, loads, and evaluate structural response.
As it concerns the first problem, outer geometry defined as a STEP model method is applied in this paper. STEP-file can support arbitrarily complex shapes, is free, easy to read, and recognized through the ISO-10303-21 standard. The outer geometry then is used to bound the lattice. The lattice itself is constructed using a standard tetrahedral mesh generator. In this study Gmsh generator [33] is used.
The output from the mesh generator consists of two matrices: connectivity matrix, later referred to as Conn.txt, and a nodal coordinate matrix, later referred to as Coord.txt. These two matrices are fed to the next step. In the next step, material and volume is associated with each strut. This enables the calculation of the response of the structure once loads are applied. The following step is to employ a multi-material physics simulator to evaluate structures' response based on the distribution of the material, loads and restraints. Standard nonlinear finite element solvers are of limited applicability in situations including large deformations, instabilities, and significant differences between elastic moduli of struts [34]. Another problem is their computational time, which makes them impractical for implementation in population-based optimization algorithms. For that purpose, VoxCAD [35][36][37][38][39][40], a multi-material physics simulator is used. Voxelyze, the underlying engine running behind VoxCAD, is a computationally efficient simulation approach based on nonlinear relaxation of structural elements. Voxelyze simulates elastic voxels based on an internal lattice of discrete points, interconnected by spring like beam elements with translational and rotational stiffness. Once it is possible to calculate deformation of the structure subjected to loads, the optimization algorithm distributes the soft and stiff lattice struts throughout the structure, trying to minimize the deviation between the actual deflection profile, and the desired deflection profile. The outline of the optimization processes are illustrated in Figure 3. To limit the search space size, the length of the cantilever beam is normalized and divided to thirteen equally spaced regions. For the problems presented in this study, this number proved to be a good balance between calculation time and ability of the structure to conform to given profiles. Finer discretization of the lattice comes at higher computational cost but would enable conforming to more complex geometries. On the other hand, a coarse discretization, with less elements, comes with increased rigidity, and thus limits the ability of the structure to conform to complex geometries. In each of these regions, the differences between the target and actual profiles are checked, and the summation of these differences defines so called root mean square error (RMSE), Equation (6), which is to be minimized. Individuals from the population with smaller RMSE values have a higher To limit the search space size, the length of the cantilever beam is normalized and divided to thirteen equally spaced regions. For the problems presented in this study, this number proved to be a good balance between calculation time and ability of the structure to conform to given profiles. Finer discretization of the lattice comes at higher computational cost but would enable conforming to more complex geometries. On the other hand, a coarse discretization, with less elements, comes with increased rigidity, and thus limits the ability of the structure to conform to complex geometries. In each of these regions, the differences between the target and actual profiles are checked, and the summation of these differences defines so called root mean square error (RMSE), Equation (6), which is to be minimized. Individuals from the population with smaller RMSE values have a higher chance of being transferred to the population of parents, and to transfer their phenotypic material to the following generations.

Algorithm Details
A population of 60 individuals is subject to optimization until either the deviation of the best individual in the population, measured by RMSE, falls below a predefined threshold, or there is no significant fitness increase in the last 50% of the evaluations. RMSE, which is the objective function, is calculated as the difference between targeted deflection profile and current deflection profile of the structure, according to Equation (6): minimize: N is for number of evaluated points along the deflection line. δ i is a deviation vector to measure the difference between targeted and current profiles, as shown in Figure 4. As indicated in previous sections, due to promotion of symmetrical load distribution through presented encoding procedures, it is enough to evaluate only one contour for RMSE minimization. In the case of symmetrical structures and loads, this suffices. For more complex distributions of loads and restraints, a minimum of two or even more profiles should be associated to the contours of the structure and included in RMSE calculation. chance of being transferred to the population of parents, and to transfer their phenotypic material to the following generations.

Algorithm Details
A population of 60 individuals is subject to optimization until either the deviation of the best individual in the population, measured by RMSE, falls below a predefined threshold, or there is no significant fitness increase in the last 50% of the evaluations. RMSE, which is the objective function, is calculated as the difference between targeted deflection profile and current deflection profile of the structure, according to Equation (6): N is for number of evaluated points along the deflection line. δ is a deviation vector to measure the difference between targeted and current profiles, as shown in Figure 4. As indicated in previous sections, due to promotion of symmetrical load distribution through presented encoding procedures, it is enough to evaluate only one contour for RMSE minimization. In the case of symmetrical structures and loads, this suffices. For more complex distributions of loads and restraints, a minimum of two or even more profiles should be associated to the contours of the structure and included in RMSE calculation. To evaluate a population and find those individuals fitter than others, a tournament selection scheme was used. Five individuals were randomly chosen from the population. The two fittest members of the five produced two children. A weight associated to the Fourier series defining material distribution was randomly selected from one and added to the other parent to create offspring. After the new population of 60 individuals was created, the mutation operator was applied.
Mutation was used to randomly change a small part of an individual and thus ensure diversity in the genetic pool and connectivity of the search space.
A randomly chosen weight of an individual was changed by 2%. In one individual, at most, 10% of weights can change through mutation.

Results
In this section, results of the proposed evolutionary algorithm for multi-material lattice structures are presented. So far, the structures are optimized to follow a predefined response, deflection after the load is applied. In all the cases illustrated in Figure 5, the left surface is clamped, the right one is free. The load was applied on the latter surface, and the bottom edge. Four cases are presented here: linear, quadratic, step, and slope discontinuity. All the deformations were multiplied ten times for the purpose of illustration. Red struts correspond to stiff material, while blue struts were for soft material. By deposition of soft and stiff struts, it was possible to evolve structures that follow given deformation behaviors with a small error. In all the results presented, the lattice consisted of 731 struts. To evaluate a population and find those individuals fitter than others, a tournament selection scheme was used. Five individuals were randomly chosen from the population. The two fittest members of the five produced two children. A weight associated to the Fourier series defining material distribution was randomly selected from one and added to the other parent to create offspring. After the new population of 60 individuals was created, the mutation operator was applied.
Mutation was used to randomly change a small part of an individual and thus ensure diversity in the genetic pool and connectivity of the search space.
A randomly chosen weight of an individual was changed by 2%. In one individual, at most, 10% of weights can change through mutation.

Results
In this section, results of the proposed evolutionary algorithm for multi-material lattice structures are presented. So far, the structures are optimized to follow a predefined response, deflection after the load is applied. In all the cases illustrated in Figure 5, the left surface is clamped, the right one is free. The load was applied on the latter surface, and the bottom edge. Four cases are presented here: linear, quadratic, step, and slope discontinuity. All the deformations were multiplied ten times for the purpose of illustration. Red struts correspond to stiff material, while blue struts were for soft material. By deposition of soft and stiff struts, it was possible to evolve structures that follow given deformation behaviors with a small error. In all the results presented, the lattice consisted of 731 struts. There were no constraints on the distribution of soft and stiff materials in the lattice, although these could be easily implemented. All the structures were evolved for five thousand generations. For the case (a) linear deformation response, the algorithm showed favorable results for this problem, being able to minimize RMSE below 10 −3 within 2000 iterations. Case (b) showed that the deflection line was constructed to have a zero slope both on the clamped and the free edge.
There were no constraints on the distribution of soft and stiff materials in the lattice, although these could be easily implemented. All the structures were evolved for five thousand generations. For the case (a) linear deformation response, the algorithm showed favorable results for this problem, being able to minimize RMSE below 10 −3 within 2000 iterations. Case (b) showed that the deflection line was constructed to have a zero slope both on the clamped and the free edge. Additionally, inflection point was enforced in the midpoint of the structure in the direction of the longest axis. To meet these criteria, we assume the deflection line is a fourth order polynomial of the form: With following boundary conditions: then the deflection line takes following form: The polynomial order is 3, which is the same as for the homogenous case, but boundary conditions are satisfied. Additionally, inflection point was enforced in the midpoint of the structure in the direction of the longest axis. To meet these criteria, we assume the deflection line is a fourth order polynomial of the form: With following boundary conditions: then the deflection line takes following form: The polynomial order is 3, which is the same as for the homogenous case, but boundary conditions are satisfied.
If we move the x location of second derivative equals zero condition away from x = 0.5, the deflection line takes the general fourth order polynomial form. The algorithm performed well on this problem, being able to reduce RMSE to approximately 5e10 −3 within 5000 iterations. Case (c) was for the step function, which proved to be hard to satisfy, especially in the region where discontinuities were located.
To minimize deflection error in this case, significant differences between elastic moduli between the two phases are required, as discussed later. Here, we have larger RMSE, approximately 2.5e10 −1 at the end, compared to the previous two cases. The last case, (d), presents the slope discontinuity, which is composed of two piecewise linear functions. The algorithm yields quality results, with the RMSE values at 10 −2 found around the discontinuity of the profile.
A typical optimization cycle of the structure is shown in Figure 6. Three phases are presented: at the beginning of the optimization, at 50% of total iterations number, and the final structure. The structure presented was optimized for negative curvature. This type of curvature is difficult to achieve, which is illustrated with the error plot. The error plot shows that difference between target and achieved is spread thorough whole structure.
If we move the x location of second derivative equals zero condition away from x = 0.5, the deflection line takes the general fourth order polynomial form. The algorithm performed well on this problem, being able to reduce RMSE to approximately 5e10 −3 within 5000 iterations. Case (c) was for the step function, which proved to be hard to satisfy, especially in the region where discontinuities were located.
To minimize deflection error in this case, significant differences between elastic moduli between the two phases are required, as discussed later. Here, we have larger RMSE, approximately 2.5e10 −1 at the end, compared to the previous two cases. The last case, (d), presents the slope discontinuity, which is composed of two piecewise linear functions. The algorithm yields quality results, with the RMSE values at 10 −2 found around the discontinuity of the profile.
A typical optimization cycle of the structure is shown in Figure 6. Three phases are presented: at the beginning of the optimization, at 50% of total iterations number, and the final structure. The structure presented was optimized for negative curvature. This type of curvature is difficult to achieve, which is illustrated with the error plot. The error plot shows that difference between target and achieved is spread thorough whole structure. The ratio between elastic moduli of stiff and soft struts is an important parameter, and it bounds the minimal RMSE the structure can achieve for the given total number of struts. This is illustrated in Figure 7. Parameter is defined as: ξ = . Figure 7. RMSE relation to different elastic moduli ratio ξ = , for negative quadratic (S1) and step discontinuity specimens (S2). The ratio between elastic moduli of stiff and soft struts is an important parameter, and it bounds the minimal RMSE the structure can achieve for the given total number of struts. This is illustrated in Figure 7. Parameter ξ is defined as: ξ = E stiff E so f t . If we move the x location of second derivative equals zero condition away from x = 0.5, the deflection line takes the general fourth order polynomial form. The algorithm performed well on this problem, being able to reduce RMSE to approximately 5e10 −3 within 5000 iterations. Case (c) was for the step function, which proved to be hard to satisfy, especially in the region where discontinuities were located.
To minimize deflection error in this case, significant differences between elastic moduli between the two phases are required, as discussed later. Here, we have larger RMSE, approximately 2.5e10 −1 at the end, compared to the previous two cases. The last case, (d), presents the slope discontinuity, which is composed of two piecewise linear functions. The algorithm yields quality results, with the RMSE values at 10 −2 found around the discontinuity of the profile.
A typical optimization cycle of the structure is shown in Figure 6. Three phases are presented: at the beginning of the optimization, at 50% of total iterations number, and the final structure. The structure presented was optimized for negative curvature. This type of curvature is difficult to achieve, which is illustrated with the error plot. The error plot shows that difference between target and achieved is spread thorough whole structure. The ratio between elastic moduli of stiff and soft struts is an important parameter, and it bounds the minimal RMSE the structure can achieve for the given total number of struts. This is illustrated in Figure 7. Parameter is defined as: ξ = . Figure 7. RMSE relation to different elastic moduli ratio ξ = , for negative quadratic (S1) and step discontinuity specimens (S2). E sti f f E so f t , for negative quadratic (S1) and step discontinuity specimens (S2).
The error plots were obtained by running the algorithm 50 times for a particular deflection profile. Two specimens were compared here: negative quadratic deflection profile (Specimen S1) and step discontinuity profile (Specimen S2).
It is obvious that step discontinuity minimizes the residual RMSE better with increase of ξ. Whereas negative quadratic profile has a sweet spot at ξ =100, and with further increase of ξ, the residual error increases. That means that there is no general rule about the nature of relation between error minimization and elastic moduli difference.
We analyzed hypothetical minimal values for ξ which will allow the desired structural deformation behaviors. We found that for simpler scenarios, i.e., linear response, and ξ > 10 to minimize RMSE to acceptable level, and ξ ≥ 1000 for most complex desired deformation profiles. We found that ξ > 10 is enough to minimize RMSE to acceptable level in case of simple deflection profiles, i.e. linear target deformation. For most complex cases, i.e. step discontinuity, a significantly larger difference in elastic moduli, ξ ≥ 1000 is required.
It is important to stress that, in theory, significant differences between elastic moduli, as high as 3000 times, are possible at the present state of the art 3D printing technology [15]. In addition to the parameter ξ, the number of lattice struts plays a role in RMSE minimization. The drawback of increasing this number is the computational cost that comes with it. The number of 731 struts used in the paper showed to be a good compromise between the quality of solutions and computational time. Fundamentally, this number does not give new qualitative insight, it is favorable to have a finer lattice compared to more coarse structures if the calculation time is not a limiting factor.

Discussion
Heterogeneity of material distribution in lattice structures was exploited to enable conforming of lattice to a set of predefined shapes. It was shown that by combining different elastic moduli, it was possible to achieve a range of predefined geometries with various complexities. Existing approaches dealing with design of functional, lattice-based materials deal with auxetics [14], optimality criteria for mass minimization [15], and design of a complex data-driven process for volume deformation control [41].
Other important findings have been made in the field of simulation tools designed for deformation analysis of lattice structures [34]. In the case of auxetic materials, careful orientation of the base units is needed to produce a handed structure with desired properties. Most of the work in auxetics deals with 2D problems. In the case of 3D structures, motions such as bending, twisting, or linear extension are considered. Recent works in 3D lattice optimization consider a grading approach to increase yield strength and stiffness of the structure through manipulation of shape and size of the struts [18].
In this study, the goal is to exploit the trait of different elastic moduli to enable conforming of the lattice structure to different 3D structures once load is applied. This approach is motivated by the current state of the art of 3D multi-material technology, which enables simultaneous printing of different material phases. The proposed population algorithm can explore the search space and evolve complex multi-material structures that deform to the predefined geometry once external forces are applied. This is not easily achieved by manually designing structures using conventional CAD techniques, or using standard topological optimization approaches, such as homogenization or shape optimization. It was pointed out that the crucial parameter for success of optimization is the ratio between the elastic moduli of stiff and soft materials, ξ.
The dependency between this parameter and residual error is neither linear nor proportionate. For some cases, residual error decreases with increase of ξ only up to a certain point, after which the structure cannot use the heterogeneity of materials as a source of further decreasing the residual error. For ξ < 10, the optimization cannot be performed, considering the problems presented in this paper and the methods used.
A procedure based on discrete cosine transform was proposed for data compression of the topology of truss structures. This increases the speed of the convergence, but also symmetry is naturally promoted in the structures. This is a favorable trait for the case of symmetrically loaded structures, and desired symmetrical deformation behaviors.

Conclusions and Future Work
The method presented in this study proved suitable for solving the problem of finding material distribution of lattices' struts with the goal of achieving desired deformation behaviors. Despite proven reliability of the used numerical simulation environment, additional analyses through both numerical and experimental validation are necessary to reliably describe the hyperelastic nature of printed structures.
The future work will include additional optimization criteria, i.e., mass minimization in addition to fulfilling deformation criteria. Another interesting problem will be to enable simultaneous shape optimization. In addition, variable lattices will be incorporated, which will affect the rigidity/compliance of the structure in certain regions, based on the structure of the lattice in that domain. The last step would be to manufacture specimens and experimentally confirm results obtained through simulations.