Validation of a nonlinear force method for large deformations in shape-morphing structures

Reducing energy and material consumption is a priority for the construction, aerospace, and automotive industries. Shape morphing addresses these concerns by broadening the band of functionality of a structure by adapting its shape to an external stimulus, such as pressure, or an internal stimulus, such as embedded actuation. This work outlines the development of an actuator placement optimization method for overdeterminate lattice structures with the objective of achieving predetermined large shape changes accurately. The deformation is modeled with both a linear and a nonlinear force method to determine their validity for large-shape change and their usefulness for the field of shape morphing. The linear and nonlinear methods are compared in four benchmark problems and two case studies relevant to the field of shape morphing. The nonlinear method is shown to achieve a level of accuracy 102\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^2$$\end{document} to 104\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^4$$\end{document} higher compared to FEM simulation, while using 23% fewer actuators and up to 77.3% less elongation of actuators, which makes it more favorable for shape-morphing applications. Two case studies for applications in aerospace and construction show that the nonlinear force method is better equipped for large shape change in overdeterminate meshed freeform target shapes and doubly curved surfaces with a high variable density. However, the nonlinear force method is less computationally efficient than the linear force method, as expected. A judicious choice of constraints can help reduce the run time.


Introduction
A global effort to reduce energy consumption and waste production is pushing engineers to rethink how to design and build structures. From the construction industry , to aerospace (Jenett et al. 2017), automotive (Daynes and Weaver 2013), and robotics (Chen et al. 2018), many have converged to a similar conceptual solution: shape morphing. A shape-morphing structure is understood to have the capacity to broaden its band of functionality by adapting its structure to an external stimulus through a shape change.
There are five main criteria to evaluate the design of a morphing structure: weight, load-carrying capacity, energy consumption, shape accuracy, and range and number of deformation modes (Wicks and Hutchinson 2004).
To tackle the criteria of weight and load-carrying capacity, much research has focused on the use of lattice structures. Their advantageous weight to stiffness ratio has motivated their implementation in airplane wings (Jenett et al. 2017), bridges , and buildings . Pellegrino (1993) and Pellegrino and Calladine (1986) showed that lattice structures are also particularly well suited for controlled deformation, as their kinematic behavior can be entirely expressed through the equilibrium matrix of the morphing structure. Shape change can be induced by inserting the appropriate number of actuators at the appropriate locations in a lattice structure (dos Santos e Lucato and Evans 2006;dos Santos e Lucato et al. 2004;Reksowardojo et al. 2019;Tian et al. 2021). Several examples of the computational approach used to place actuators in a structure for shape morphing are given in Table 1.
They address some of the remaining criteria of shape morphing, but not entirely. Methods that assume small deformations (strain ≤ 1%) significantly impair their accuracy for large shape change, or reduce the range and number of deformation modes. The method presented by Reksowardojo et al. (2020) addresses all criteria the most completely but limits its application to five lattice cells and 26 variables (see Table 1). The use of a linear or nonlinear method impacts all the criteria of shape morphing but have not yet been compared directly.
The first contribution of this work is the strain-based nonlinear force method, which models how a strain input affects a structure as a whole in the optimization of the placement of linear actuators in a structure to achieve a target shape. The method is adapted from Reksowardojo et al. (2020), whose primary contribution is to counteract the effect of external loads through shape change. The second contribution is the inclusion of the optimization of the number and placement of actuators in both methods. The third contribution is the quantitative comparison of both methods to show that the proposed nonlinear force method is better suited for determining the necessary actuators and their movement for a structure to morph into several target shapes. Most work in the literature has focused on using linear FEM because it requires less computational effort and can be implemented directly in an optimization method. However, it limits the deformation range of a structure and loses accuracy in structures with a large number of states of self-stress.
Both methods presented are benchmarked for four main deformation modes and a similar number of cells and variables against the examples given in Table 1 with respect to the shape-morphing criteria: energy consumption (or number of actuators required), shape accuracy, and range and number of deformation modes achieved.
Finally, both methods are applied to two case studies ubiquitous in the field of shape morphing: a 2D NACA wing profile and a 3D morphing panel. These case studies demonstrate the contributions of this work by extending the field of application to overdeterminate meshed freeform shapes and doubly curved surfaces with a high variable density. The structure of this paper is as follows: first, the theory behind the linear and nonlinear force methods to determine the displacements resulting from the inclusion of actuators is explained in detail. The implementation of the actuator placement optimization using each method is detailed. Then, the benchmark problems, the case studies, and their evaluation parameters are introduced. The benchmark deformations are compared to FEM simulations conducted with Abaqus CAE. The linear and nonlinear force methods are evaluated first with respect to the benchmark problems and then with respect to the case studies. Finally, the results are discussed first in the context of this work and then in the greater context of shape morphing.

Actuator placement optimization
The search for the appropriate number and location of actuators has been addressed using several numerical approaches.  (2021) use a statically and kinematically determinate Kagome lattice structure to produce a hinging and twisting plate. Static and kinematic determinacy allow to model the effect of a set of actuator elongations as the linear combination of single-actuator elongations. This method is extremely quick once the effect of each actuator has been registered in a matrix, but it is limited to small deformations and determinate structures. Senatore et al. (2019) use a combination of force and displacement sensitivity matrices and load path redirection to first determine the most efficient actuator locations, and then the internal force method to optimize their lengths. This method assumes small deformations: its cases studies achieve a maximum element elongation of 0.34%. Reksowardojo et al. (2020) use simulated annealing (SA) and the Tanimoto index, a shape similarity criterion, to determine actuator locations. The actuator elongations are determined using force and displacement sensitivity matrices and optimized with a Newton-Raphson scheme. Force and displacement commands are determined using a nonlinear method, which allows for large element elongations, 67.2% at the most. Wang et al. (2020) and Tian et al. (2021) establish control sequences based on a nonlinear kinematic method using constant curvature, segmentation method, and beam theory. They achieve an exact curvature change in a 2D morphing panel in three steps using two lattice cells by elongating the outer edges using linear slide modules, thus severely, restricting their achievable design space.

Modeling of strain-induced deformation using the FM
Most work in structural analysis focuses on the displacement based finite-element method (FEM) (Luo and Lu 2006). However, the works of Pellegrino and Calladine (1986), Pellegrino (1993), and Patnaik (1973) are the basis for formulations of the force method (FM) that have inherent advantages in modeling adaptive lattice structures. The main building blocks for the FM are the equilibrium equations and the compatibility equations. The equilibrium matrix A and its subspaces explicitly link the structural and kinematic properties of a lattice structure. Actuator elongations can be formulated as imposed strain distributions so internal forces and nodal displacements can be calculated directly. However, the FM requires a first empirical step of eliminating redundant lattice members and replacing them with external loads. Also, the equilibrium matrix is rectangular and potentially asymmetric for overdeterminate systems and is, thus, not sufficient to solve for internal forces directly. Pellegrino (1993) and Patnaik (1973) developed two formulations to avoid these issues. The approach of Pellegrino (1993), building on earlier work by Pellegrino and Calladine (1986), is based on the singular value decomposition (SVD) of A and the interpretation of its subspaces. Starting with the equilibrium equations and the compatibility equations, he formulates a method to derive static and kinematic determinacy of a structure, s and m. Subdividing the SVD of A based on s and m, he derives the nodal displacements and element forces resulting from imposed strain distributions and external load. For a statically and kinematically determinate structure ( s = 0 and m = 0 ), the solution is direct and exact. This formulation is only valid for small linear displacements and is extended by Luo and Lu (2006) to the nonlinear force method (NFM) by using a Newton-Raphson scheme. This formulation of the FM has recently been used by Reksowardojo et al. (2020) and Kim et al. (2016), focusing on keeping structures within serviceability limits, compensating for vibrations, or detecting damages.
The approach of Patnaik (1973), the internal force method (IFM), avoids the expensive computation time of SVD by compensating for the rank deficiency of A. It builds a new matrix = [ red | ⋅ ] T , where the equilibrium matrix is reduced to free degrees of freedom (DoF) red and is completed with the multiplication of the compatibility matrix C with the local flexibility matrix G. The additional equations ensure that the desired element elongation is compatible with an imposed strain distribution.
The concept of sensitivity analysis in adaptive lattice structures is introduced in Senatore et al. (2019), Senatore and Reksowardojo (2020), and Yuan et al. (2016). The force and displacement actuation sensitivity matrices F and u contain the element forces and nodal displacements triggered by a unitary length change for each of the structure's elements in turn. Element forces and nodal displacements caused by a given element elongation vector x are given directly by F ⋅ and u ⋅ for small displacements. F and u can be calculated using both the FM and the IFM, but their calculation is a lot more efficient with the FM formulation. The importance of the sensitivity matrices in the optimization of actuator placement in lattice structures motivates the use of the FM formulation of Pellegrino (1993) and the NFM formulation of Luo and Lu (2006) in this work.

Method
This section describes the theory and implementation for the actuator placement optimization assuming first that the effect of elongating elements is linear and then that it is nonlinear. Both the linear and nonlinear methods are based on previous work by Pellegrino and Calladine (1986) and Pellegrino (1993) and their definition of kinematic and static determinacy.

The equilibrium matrix
The equilibrium matrix A of a lattice structure contains all information pertaining to its static and kinematic behavior. It relates external loads to internal loads directly. The SVD of this matrix separates the information into three matrices: where A ∈ ℝ n D ×n E , U ∈ ℝ n D ×n D , W ∈ ℝ n E ×n E , and ∈ ℝ n E ×n E . r m , r s , and r are the left singular vectors, right singular vectors, and singular values of matrix A, respectively. contains information on the kinematic behavior of , and on the static behavior. r ∈ ℝ n D ×(n D −m) contains the load components in equilibrium with the stress components in r ∈ ℝ n E ×(n E −s) multiplied with the corresponding singular value in , so it is the basis of the row space of , R( ) . Since the components of r are in equilibrium, they do not cause any displacements. s is the basis of the null space of , N( ) . m ∈ ℝ n D ×m contains the load components that are not in equilibrium and, thus, cause displacements and form the basis of the inextentional mechanisms of the structure. Similarly, s ∈ ℝ n E ×s contains the s linearly independent states of self-stress of . s can also be interpreted as the incompatible strain basis.
A structure is statically determinate if it contains no states of self-stress, or s = 0 , and is kinematically determinate if it contains no mechanisms, or m = 0 . The values of m and s are calculated directly from A: and where n E is the number of elements in the structure, or the number of rows of A; n D = dim ⋅ n N − n k is the number of nodes multiplied by the dimension of the structure minus the number of kinematic constraints. r is the rank of A. If is kinematically determinate, m is empty. If is statically determinate, s is empty.

The force method for imposed strain
The force method used in this work is based on the work of Pellegrino (1993), Senatore et al. (2019), Senatore and Reksowardojo (2020) and is referred to as FM. The FM is based on two sets of equations: the equilibrium conditions and the compatibility conditions. The force equilibrium conditions at a structure's nodes are where is the internal load vector and is the external load vector. The compatibility conditions as defined for the SVD method by Pellegrino (1993) ensure that compatible strain should vanish through orthogonality to the incompatible strain basis s : or where ∈ ℝ n E ×n E is the element flexibility matrix and Δ ∈ ℝ n E ×1 is the element elongation vector. Through the principle of Virtual Work, the compatibility conditions can be rewritten for statically and kinematically determinate systems as Pellegrino and Calladine (1986): so where are nodal displacements. The FM, thus, allows the calculation of element elongations directly from the nodal displacements of a target deformation for kinematically and statically determinate structures, assuming that the effects of the elongations are linear.

Sensitivity matrices
For indeterminate structures and assuming small deformations, element elongations can be calculated through the force and displacement sensitivity matrices F ∈ ℝ n E ×n E and u ∈ ℝ n D ×n E . F relates the structure's flexibility and states of self-stress to determine reactive forces. u relates the structure's flexibility, force sensitivity matrix, and the pseudoinverse equilibrium matrix to determine displacements ( Senatore et al. (2019)). F i and u i contain the effect of the unitary length change of the ith element of Δ , excluding the external load . Another advantage of the FM formulation is that the two matrices can be calculated directly rather than with a loop (Yuan et al. 2016;Senatore and Reksowardojo 2020):

Optimization objectives
The objective of this work is to find the actuator placement and elongations in a lattice structure that best approximate a given target shape, assuming that all elements can be actuators. For indeterminate structures, this problem translates into the following optimization: where and are the maximum contraction and elongation for each element, and is a vector of length n E × 1 . x i is the value of the elongation of the ith element, so x i = 0 if the element is not an actuator. The quality of a solution x is evaluated with the Tanimoto index, a shape similarity criterion used in machine learning (Tanimoto 1958;Reksowardojo et al. 2020): where Q i is the Tanimoto index for the unitary elongation of element i, c is the controlled deformation, and target is the target deformation. The closer the Q is to one, the better the shape approximation. In general, the lowest objective function value corresponds to the highest Tanimoto index value but not always since the latter is both a shape and a node similarity criterion, as opposed to the objective function which is only a node position similarity criterion.

Linear method
The linear method assumes small deformations of the elements and large nodal displacements. Elements are allowed to elongate and contract up to 100% of their initial length to compare the validity of the model for large-scale shape morphing to the nonlinear method. Figure 1a shows the optimization framework used for the linear method: after calculating the sensitivity matrices as explained above, the actuator placement is optimized using element displacement efficacy assuming that the number of actuators is predetermined, and then the actuator lengths are optimized to achieve the desired target shape target .

Actuator placement
The placement of the actuators within the structure is determined using the concept of placement efficacy from Senatore et al. (2019). The efficacy is a measure of the contribution of each individual element to the desired displacement. The closer the efficacy of an element is to 1, the more it contributes to the desired displacement. It is evaluated as follows: u|cDOFs j is built by extracting the rows from u corresponding to the jth controlled DOF. Δ i is the unitary elongation vector of the ith element. cDOFs target is the target vector reduced to the controlled DOFs.
The necessary number of actuators is determined as the number of actuators needed to solve Eq. 14 exactly. For the linear method, the actuator positions are chosen as the n act elements with the highest efficacy.

Actuator length optimization
Assuming small displacements and no external loads, Senatore and  show that element elongations to achieve a target shape can be directly calculated with where the + indicates the Moore-Penrose inverse. However, the quality of the results depends on how close r and u are to singularity. If the norm of the difference between target Page 6 of 17 and 0 , the nodal displacement caused by , is smaller than the convergence criterion e Lin , then x is already the optimal solution, so the number of objective function calls is 0 and the Tanimoto index Q = 1 . If the norm of the difference Fig. 1 Optimization method for the linear method (left) and nonlinear method (right); the NLE solver is described in the Replication of results section between target and 0 is larger than e Lin , the element elongations are optimized using Eq. 11 reduced to the controlled DOFs and the known actuator positions: The starting point 0 is previously calculated with Eq. 14.
The quality of opt is determined by the resulting Tanimoto index.

Nonlinear method
The nonlinear method presented in this work can be used for large deformations ( >> 1% ). The nodal displacements caused by element elongations are computed using a nonlinear solver based on the work of Pellegrino (1990), Luo and Lu (2006). maps the initial structure and given elongations to the resulting nodal displacements, element forces, and Tanimoto index Q: where and are the cross-section vector and elastic modulus vector for all elements. Elements can elongate and contract up to 100% of their initial length to assess applicability of the method to large-shape changes. Figure 1b shows the optimization method used for the nonlinear method introduced in this work. A shape-morphing lattice structure should contain as few actuators as possible to not compromise its structural soundness. Reksowardojo et al. (2020) define that The first optimization shown in black in Figure 1b starts at n act 0 = n cDOFs + s . It is implemented to calculate the necessary number of actuators, first using simulated annealing (SA) to determine the actuator positions given unitary elongations, and using IPOPT to optimize the actuator elongations given the positions. SA is chosen due to the combinatorial nature of the actuator placement problem. In the first step ( i = 1 ), the placement and lengths are optimized r times. In the second step ( i = 2 ), the best result from i = 1 Q best is saved as Q old , and actuators are added so n act = n act + n E 100 .
⎧ ⎪ ⎨ ⎪ ⎩ u = target for n act = n E u ≃ target for n cDOFs + s < n act < n E u ∼ target for n act < n cDOFs + s . The placement and lengths are optimized again for the new n act . Convergence is achieved if the new best solution Q best is worse than Q old , or Q best Q old < 1 . The first optimization method is shown in Algorithm 1.
The second optimization method calculates the best actuator placement and elongations, given n act opt . It loops through the same SA and IPOPT algorithms as the first optimization and converges either when the Tanimoto index Q best exceeds the convergence criterion Q conv or when the maximum number of iterations is achieved. For each step, the Tanimoto index, the value of the objective function, and the iteration number are stored in a matrix if they exceed the previous best result.
This work modifies the method reported by Reksowardojo et al. (2020) by optimizing an imposed internal strain to achieve a target shape deformation rather than minimizing any deformations caused by an external load p. It further enhances the method by introducing loops within both optimizations. This meta-heuristic approach enables both optimization methods to avoid local optima in this nonlinear problem.
The linear method is used to determine the nodal displacements given the element elongations x within the objective function of the IPOPT step, given in Eq. 15, for computational efficiency. The results are then adjusted using a nonlinear method before calculating the Tanimoto index and checking for convergence. The details of the implementation of the actuator placement using SA and the elongation optimization using IPOPT are given below.

Actuator placement using SA
The nonlinear nature of the system nullifies the assumption that the n act most efficient elements constitute the best actuator placement solution, since it takes into account only the efficacy of a single actuator at a time rather than actuator combinations. Instead, the efficacy is evaluated using the Tanimoto index for the unitary elongation of each element. The number of possible actuator combinations increases exponentially with the number of elements, so an exhaustive search is impractical. Since the Tanimoto index still only describes the efficacy of single elements but not their combination, the vector Q containing the Tanimoto indices for each element elongation is assumed to represent a probability distribution for elements to contribute to target .
The distribution and its pseudo-inverse, pinv , are used in a SA stochastic search to determine the best actuator placement given an initial n act and assuming unitary elongations. The problem size is reduced by setting the probability distributions to zero for elements that link two pinned nodes, as they cannot elongate. The method, which mimics the cooling schedule of molten metals through annealing, was first formulated by Kirkpatrick et al. (1983). In the beginning, as the temperature T is high, the algorithm is likely to accept inferior neighboring solutions and avoids local optima. As T reduces, the algorithm is likely to reject inferior neighboring solutions, and the search is increasingly focused towards an optimum.
The goal is to maximize the Tanimoto index for a combination of actuators. The objective function is where is a penalty. The starting point 0 is chosen as the n act elements with the highest Q j . Neighboring solutions j are generated by drawing without replacement from the n E − n act remaining elements using the probability distribution and replacing certain actuator elements using the pseudo-inverse distribution pinv . In each iteration, the number of elements replaced in j is chosen at random, and the Tanimoto index of the combined unitary elongation of the new actuator set is evaluated. The acceptance probability A Q for a new solution is dependent on T: is draw at random from a uniform distribution U(0, 1) . The acceptance probability A for a new penalty is also dependent on T: As T sinks, so does the acceptance probability of a lower penalty value, and so the search is increasingly focused.
T is updated n temp times, and the search is repeated n rep times for each temperature. T k , the temperature at step k, is defined as follows: The best solution found so far is stored. Convergence is achieved when the solution no longer improves after several iterations. The resulting solution is the actuator placement that scored the highest Tanimoto index. Further implementation details on this method can be found in Wah and Wang (1999) and Reksowardojo et al. (2020) and in the Replication of results section of this work.

Actuator elongation optimization using IPOPT
Once the actuator placement has been determined assuming unitary element elongations using SA, the element elongations are optimized using IPOPT. Since the objective function optimization given in Eq. 15 is estimated using the linear model for computational efficiency, the process is iterative. First, Eq. 15 is solved using IPOPT. Then, the objective function value and Tanimoto index are calculated for the resulting element elongation j using the nonlinear model. Convergence is achieved either if the objective function value is smaller than e NL or if the Tanimoto index is lower than in the previous step. If convergence is not achieved, the initial actuator elongations of the next step are set to j and the coordinates, equilibrium matrix, and sensitivity matrices are updated to reflect the element elongations of j . The boundary conditions are also updated so the total element elongation boundary remains the same. The element elongation optimization is shown in detail in Algorithm 2. .

Benchmark problems
Both methods are tested on four 2D target deformations given in Table 2 that are benchmarked against results using Abaqus CAE 6.14: 1. linear deformation, 2. concave up, 3. concave down, and 4. sine, the transition between the two. The method parameters used in the optimization method as well as the simulation are given in Table 3. All problems have static determinacy s = 1 and kinematic determinacy m = 0 . The quality of both methods is assessed using the Tanimoto index of the optimized results, the average error of their nodal displacements with respect to the simulation, the number of actuators needed for each target deformation, the absolute largest strain, the number of iterations, and the number of objective function evaluations.

Case studies
The methods are tested on two case studies: in 3D, for complex doubly curved deformations of a tetragonal lattice structure plate, and in 2D, on a NACA 2640 airfoil, using an Abaqus CAE generated tetragonal mesh.

Implementation details
The work presented above is implemented in MATLAB R2020a and verified with Abaqus CAE 6.14-1 using a static general step and NLGEOM. This problem is solved using an interior point optimization (IPOPT) made available by (Wachter and Biegler, 2004) with default settings and maximum iteration number set to 250. The optimization is run on the ETH Euler cluster using ten cores with 64-core AMD EPYC 7H12 processors.

Results
This section first details the results for the benchmark problems of both methods and then the results from the case studies.
For the linear method, the convergence criterion e Lin is set to 10 −3 and the maximum number of iterations is 100.
For the nonlinear method, the convergence criterion Q conv is set to 0.999 and the maximum number of iterations is 100. Preliminary studies of the stored results matrix showed that setting the maximum number of iterations higher did not lead to a better result for the structures analyzed in this work. The number of loops r in the actuator number optimization is set to five. Preliminary studies showed that increasing the number of loops to ten or twenty only lead to better results in less than 10% of the cases; the number of loops is, thus, set to five for computational efficiency.
Preliminary studies where results were calculated exhaustively for n cDOFs + s ≤ n act ≤ n E showed that, for problems with values of n E > 100 , the ratio Q best Q old increased every n E 100 steps. The size of a problem with n E 100 also motivated the reduction of the number of steps. The convergence criterion for the actuator length optimization e NL is set to 10 −6 .

Benchmark
The benchmark results for the two methods are detailed in Fig. 2 and Table 4. Comparing the LE and the NLE method, the Tanimoto indices of all results vary by less than 0.2% from the target shape, but the error with respect to the simulation are 10 2 to 10 4 orders of magnitude larger for the LE. This is also visible when comparing the left and the right columns of Fig. 2a and b, especially in the fourth benchmark problem. The LE method does not take into account internal stresses caused by other element elongations, which in fact cause the structure to warp significantly. Further tests show that the LE method can only achieve similar levels of accuracy if limited to less than 10% element elongation, which in turns limits the achievable macroscopic deformations and drastically increases the number of actuators. The NLE method needs on average 23% less actuators than the LE method. Its actuators elongate between 1.8% and 77.3% less. The actuator placements are analogous, but better calibrated for the NLE method. Comparing for example the sine deformation in Fig. 2.4, the contracting and elongating elements are placed similarly for the Fig. 2 Benchmark results for both methods (a and b) and four deformations (1 to 4) with static determinacy s = 1 ; for each deformation type: upper left-actuator placement, blue is elongation, red is contraction; lower left-deformed structure from the method, the red circles indicate the target positions for moving nodes, the black circles indicate fixed nodes; right: axial stress for the deformed structure from the simulation; the stress scales on the right of each column apply for all simulation results concave down and up portions of the structure, but the NLE method makes more efficient use of the placement, using two rather than eight actuators in both portions, which can lead to a more structurally sound and energy efficient structure. The LE method has the advantage of converging faster, needing roughly 200 times less objective function evaluations and on average 4.64 times less run time. The run time for the LE is under one minute whereas the NLE method needs up to three minutes. The average run time per objective function evaluation is 7 × 10 −2 s for the LE method and 2.3 × 10 −2 s for the NLE method, so the NLE method actually iterates faster. The convergence diagrams for the nonlinear actuator number optimization of the benchmark problems are shown in Fig. 3, where the exact results given in Table 4 are indicated in red. The aim of this work is to obtain the best compromise of shape similarity to number of actuators. As shown in the figure, adding actuators beyond the convergence criteria given in Sect. 3.6 and at the beginning of this section offer little to no improvement.

3D panel
The results for the 3D panel for the two methods are detailed in Fig. 4 and in Table 5.
The NLE doubly curved deformation requires nearly four times as many actuators as the LE method. The actuator pattern in Fig. 4a shows that the LE method only deforms locally to accommodate the target shape, whereas the NLE method presents a wider actuation response. The pattern of contracting and elongating elements is, however, apparent in both cases.
The corner lift deformation differs from the doubly curved deformation significantly. For the first and only time, the maximal strain of the actuators for the LE method is significantly lower than for the NLE method. However, it requires nearly three times as many actuators overall, so the lower strain is compensated by a full body response. The NLE method requires significantly more objective function calls and on average 6.38 times more run time, but its resulting Tanimoto index is perfect. Setting lower convergence limits could cut down the number of function evaluations, although it would affect the accuracy. The highest run

NACA airfoil
The results of the meshed NACA airfoil for the two methods are detailed in Fig. 5 and in Table 6. This case study has five times more nodes than the benchmark problems and twenty-five times more elements. The increase in the elements to nodes ratio considerably increases the numbers of states of self-stress in the structure which affects the results. As opposed to the benchmark problems, the LE method achieves a Tanimoto index of only 0.793, compared to a nearly perfect 0.993 for the NLE method. This is clearly visible in Fig. 5a for both deformations, where the method results overshoot the desired deformation. However, the increase in shape accuracy comes with a cost, since the NLE method needs 41.7% and 5.5% more actuators for

Discussion
This section first reviews the results presented and shows the superiority of the NLE method quantitatively. It then compares qualitatively to examples of shape-morphing lattice structures from literature and discusses the relevance of this work in the broader field of shape morphing. The results from the benchmark problems clearly show the advantages of the NLE method. With respect to the benchmark problems, the solutions for the NLE method are more accurate and need on average less actuators, which is an advantage for a structure's soundness. The doubly curved 3D case study is the only instance where the NLE method needs more actuators at a slightly higher maximum strain. Using the actuator location combinations efficacy calculated with the SA rather than single-actuator efficacy combined with the NLE method consistently achieves better results with less actuators and less elongation. The only benchmark problem where the NLE method needs more actuators than the LE method is the sine deformation, which is also where the LE method deviates the most from simulation. However, the accuracy comes at the cost of computational efficiency, with calculation times roughly four times longer than the LE method. The LE method does not take into account internal stresses caused by other element elongations, which in fact cause the structure to warp significantly. Further tests show that the LE method can only achieve similar levels of accuracy if limited to less than 10% element elongation, which in turn limits the achievable macroscopic deformations and drastically increases the number of actuators.
These differences are made more obvious with an increase in number of variables and structural complexity. Adding elements, and thus, states of self-stress, impacts the computational effort: multiplying the number of variables by 25 results in a calculation time roughly 100 times longer.
The effect on the run time per iteration is different for the LE and NLE optimizations. In the NLE optimization, increasing the degree of overdeterminacy from 1 to 207 results in a multiplication of the run time per iteration by a factor of 1,374, but a further increase to s = 509 in the 3D panel case study only increases the run time per iteration by a factor of 3.60. In the LE optimization, the same increases correspond to a factor of 5.14 for s = 1 to s = 207 and a factor of 154 for s = 207 to s = 509 . So, although the LE method is generally more efficient, it loses its advantage with an increasing degree of overdeterminacy. The NLE methods is thus more appropriate for very large complex structures.
When comparing the results from the two case studies, the 3D panel deformation requires up to 24 times less calculation effort than the NACA profile for the LE method, whereas it requires 2.5 times more calculation effort for the NLE method. The LE method, where the impact of the states of self-stress is less significant, is impacted exponentially by the number of elements and the deformation complexity. However, the NLE method needs more objective function evaluations in the NACA case where there are less variables, but a higher variables-to-states of self-stress ratio. In the NACA case study, the constraints are often on both ends of elements, which creates inherent self-stress, and further increases the number of objective function evaluations needed to solve for the target deformation, which shows that the position of constraints impacts the efficiency of the methods greatly.
However, the case studies differ from the benchmark in the number of actuators needed for both methods. Since the LE method does not take into account the internal stresses caused by element elongations, it underestimates how many actuators it will need to achieve the target deformation, which in turn affects the Tanimoto index. This problem is particularly apparent in the NACA case study, where several elements are self-stressed by the constraints. The NLE method deforms incrementally and includes the induced states of self-stress, which increases the number of necessary actuators significantly while conserving the accuracy of the model. On the other hand, the LE method overestimates the number of actuators needed for the panel "corner lift" case by a factor of four, and delivers a structure that would be significantly harder to build that the NLE method. The case studies also differ from the benchmark in that they need 10 to 100 times fewer objective function evaluations. This is explained by two factors. First, the heuristic for Optimization 1 takes steps of n act = n E 100 , which is 1 for the benchmark, 15 for the 3D panel, and 5 for the NACA profile. This means that, for example, the Corner lift deformation of the 3D panel only iterates twice before n act converges. Second, the ratio of deformation to structure size is much smaller in the case studies, so they require a smaller actuator density than the benchmark to achieve the target shapes.
The NLE method supplemented with the stochastic actuator number optimization (optimization 1) and length optimization (optimization 2) addresses the five main criteria of shape morphing better than the LE method: weight, loadcarrying capacity, energy consumption, shape accuracy, and range and number of deformation mode. As in comparable work given in Table 1, weight and load-carrying capacity are addressed in the use of lattice structures. By using actuators more efficiently both in number and maximal strain, the NLE method shows better load-carrying capacities and better energy consumption. The shape accuracy and deformation range are similar for both in simple problems with low levels of overdeterminacy, but the NLE method is better equipped for a range of complex target shapes and high levels of overdeterminacy.
Comparing back to the examples given in the Introduction, the NLE method addresses all criteria where many addressed only one or two. Comparing further to Reksowardojo et al. (2020), on whose work some of the method is based, the NLE method is shown to achieve a larger shape range and to work with a more complex design space.
The extension of the design space and target shape range to highly overdeterminate freeform meshes and variable dense 3D problems of the case studies show the promise of the optimization method proposed in this work for the field of shape morphing. The nonlinear method can be applied to problems with any arbitrary lattice layout, variable number, and target shape. However, an increase in a problem's complexity is mirrored by an increase in its computational effort. Further work should focus on finding an acceptable compromise between accuracy and computational effort with certain simplifications of the NLE method.

Conclusion
This work describes a linear and a nonlinear method to estimate the shape-morphing capacities of an overdeterminate lattice structure by optimizing the number and position of linear actuators within it. It compares the methods quantitatively in four benchmark problems and two case studies relevant to shape morphing. It shows that the nonlinear method is the most appropriate for large deformations of complex structures and leads to results of up to 10 4 higher accuracy than the linear method. By implementing a new stochastic optimization approach to minimize the number of actuators in the nonlinear method, it delivers accurate results that are more structurally sound than the results from the linear method. The case studies validate the method for highly overdeterminate, meshed freeform target shapes and doubly curved surfaces with a high number of variables, and indicate the importance of the choice of constraints in these systems. The nonlinear method can, thus, be applied to problem sizes with a number of elements on the order of 10 3 and even larger, target shapes with varying curvature, and high degrees of overdeterminacy. However, the efficiency of the nonlinear method is affected by the number of variables, the constraints, and the degree of overdeterminacy. The case studies show the potential for the method to optimize the design of complex, active structures applied in aerospace and construction. Future work should further test the limits of this method in terms of number of variables versus achievable deformation, and the results should be verified experimentally.

Nonlinear method implementation
The NLE method employed in this work is based on Luo and Lu (2006). This appendix describes its implementation in Matlab R2020a in detail.
The process iterates in order to reduce the residual strain and the residual force, both of which are possible inputs. So, given the strain vector and/or the external force vector F, the process outputs the elements strains and nodal displacements. In this work, the equilibrium matrix is extended to include support reactions, since element elongations affect the supports directly. The formulation of the extended equilibrium matrix ext is detailed in Senatore and Reksowardojo (2020). The matrix is built as follows: where el is the usual equilibrium matrix, as defined by Pellegrino and Calladine (1986), and sr contains the effect of the support reactions. If the support reactions are aligned with the main axes of the system, as they usually are, this matrix contains all zeros except for the constrained nodes, where they contain ones for each constrained direction. The flexibility matrix G is also necessary for the NLE method, and so also needs to be extended for the support reactions. Since the supports are considered to be infinitely stiff elements, ext contains n fDOFs additional rows and columns, all set to zero.

Simulated annealing implementation
In this work, a simulated annealing (SA) algorithm is implemented to optimize the actuator placement given n act . Its main characteristics are given in Sect. 3.6.1, and the details of its implementation are given below. n temp and n rep are the number of temperature steps and the number of iterations at every step respectively. They are set to 30 and to n E . The initial actuator layout 0 is determined using the LE method. The initial penalty is chosen at random from U(0, 1) . The objective function is defined in Eq. 19 and is defined here as J.
In this work, the convergence criterion is based on the variability of the entries of . If they have to changed over four temperature increments, the algorithm outputs J best . otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.