An inverse method for optimizing elastic properties considering multiple loading conditions and displacement criteria

Signiﬁcant research effort has been devoted to topology optimization of two- and three-dimensional structural elements subject to various design and loading criteria. While the ﬁeld of topology optimization has been tremendously successful over the years, literature focusing on the optimization of spatially-varying elastic material properties in structures subject to multiple loading states is scarce. In this article, we con-tribute to the state of the art in material optimization by proposing a numerical regime for optimizing the distribution of the elastic modulus in structural elements subject to multiple loading conditions and design displacement criteria. Such displacement criteria (target displacement ﬁelds prescribed by the designer) may result from factors related to structural codes, occupant comfort, proximity of adjacent structures, etc. In this work, we utilize an inverse problem based framework for optimizing the elastic modulus distribution considering N target displacements and imposed forces. This approach is formulated in a straight-forward manner such that it may be applied in a broad suite of design problems with unique geometries, loading conditions, and displacement criteria. To test the approach, a suite of optimization problems are solved to demonstrate solutions considering N = 2 for different geometries and boundary conditions.


Introduction
Determination of optimal structural geometries, shapes, connectivity, cross-sectional area, and void locations are central themes in topology optimization. In cases where additional degrees of freedom for structural optimization are desirable, one may additionally con-sider the optimization of macroscopic material properties (such as elastic properties), which is the central theme of this article. With this said, the ability to optimize structural topology -alone -affords designers substantial freedom by reducing structural weight, increasing adherence to architectural constraints, and providing an additional means to meet numerous other design criteria [1][2][3]. Significant research has been devoted to topology optimization over the last thirty years, often aiming to optimize geometrical aspects of structural elements (voids, connectivity, thickness, etc.). In this effort, various topology optimization approaches including the well-known "solid isotropic material with penalization" (SIMP) regime [4][5][6], level-set methods [7][8][9][10][11], phase-field methods [12][13][14], and other emerging regimes have had much success.
Complimenting topology optimization is the optimization of elastic material properties (herein referred to as simply material optimization), wherein the distribution of elastic material properties is optimized with respect to the design criteria [15]. Material optimization may be broadly viewed as a part of "free material design" (FMD), which often aims to maximize structural stiffness by optimizing the material and material distribution given some design criteria and constraints [16]. Classically, these criterion are incorporated by prescribing the cost of design; generally, the cost Π is is defined by some integral over the domain of interest Ω as follows where Φ is a scalar function and C is an elastic tensor in the Cartesian coordinate system x. Arguments incorporated with the scalar cost function Φ(·) vary upon the design criteria, optimization variables, and prescribed constraints. In cases where only the material thickness, voids, topology, etc. are optimized, the values of the elastic constants are generally distributed equally within Ω.
The attractiveness of material optimization is that it affords designers additional flexibility by increasing the degrees of freedom within Ω -since the models permit inhomogeneity in the elastic parameters.
Recently, there has been increased interest in material optimization (e.g. [17][18][19][20][21]); this may, in part, be inspired by emerging engineered materials such as functionally-graded composites, 3D printed materials, fiber-reinforced cement-based composites, and material produced using additive fabrication [15,22]. One potential advantage of these engineered materials is that their material properties, such as the elastic modulus, may be inhomogeneous across the structural geometry. Therefore, the realization of such technologies may promote the usage of material optimization in a large suite of structural applications.
Computational material optimization regimes, often inter-weaved within FMD frameworks, have been well-studied since their origins in seminal works (e.g. [23,24]), where the authors developed numerical approaches for optimizing materials to achieve minimum compliance. Soon thereafter, researchers utilized regimes employing material and stress constraints (i.e. interior-point methods and penalty-barrier methods), demonstrating improved efficiency in their optimization regimes [25][26][27][28]. One aim of constraining the stress during the material optimization process is to prevent structural failure occurring at high stress levels. However, later works showed the difficulties of using stress constraints [29,30]. The difficulties in utilizing stress constraints often resulted from defining suitable failure criterion and technical complexities in programming stress constraints nested in optimization algorithms [31]. In overcoming these challenges, authors [31][32][33] developed augmented Lagrangian methods, demonstrating their effectiveness in overcoming obstacles resulting from constraints. Extending the extensive body of work related to constrained materials optimization, researchers in [34,35] advanced and broadened the state of the art by proposing generic constraints (including displacement and local stress constraints) dependent on design and state variables. Researchers in, e.g. [36][37][38], have also re-cently demonstrated the effectiveness of topology optimization regimes incorporating stress constraints coupled with shape sensitivity functions, level-set methods, geometric projection methods, and much more.
Within recent years, substantial research in material optimization has been conducted. Some of the major contributions during this time include publications in isotropic material design (IMD) [39], cubic material design (CMD) [40][41][42], and Young's modulus design (YMD) [15]. Broadly, the use of these material optimization methods, among the many others mentioned, may be directly applied in conjunction with manufacturing processes utilizing, for example, programming of cellular materials [43], truss-like cellular structures [44,45], cement-based materials [46,47], and functionally-graded materials [48,49]. Moreover, much attention has also been given to optimization problems where the primary focus is not purely related to elasticity. Some examples include the optimization of transport properties in porous materials [50][51][52], optimization of vessel geometry to promote fluid flow [53][54][55][56], and topology optimization for heat conduction problems [57][58][59].
Although the body of work in material/topology optimization is quite substantial, often many approaches consider only a single load case. In practice, few structural elements experience a unique loading condition through their life cycle. This realization was formally addressed in [16], where the authors proved the existence of a solution to the multi-load optimization problem and a regime for solving the problem with respect to the distribution of the materials as well as the material properties. Since this seminal article, works aiming to optimize structures under multiload conditions have been noted (e.g. [60][61][62][63][64][65]). In general, however, literature focusing on the optimization of spatially-varying elastic material properties in structures subject to multiple loading states is scarce.
In this work, we address this gap in research by developing a generic and straight-forward computational framework aiming to determine an optimal distribution of the elastic modulus considering multiple loading conditions and displacement criteria (target displacement fields prescribed by the designer). We reinforce that the aim of the proposed approach is not to optimize the topology or geometry of the structure of interest; rather, we focus solely on the distribution of the elastic modulus. Concisely, in this article we aim to accomplish the following: 1. Develop a straight-forward and generallyapplicable inverse problem based framework for optimizing the elastic modulus distribution considering multiple loading criteria and imposed forces. 2. Determine the effects optimization parameters, the geometry of the structural members, and imposed forces/displacement criteria have on optimized solutions. 3. Analyze the effects that discretization and imposed constraints have on the regime's performance and optimized solutions.
The article is organized as follows. We begin by presenting the mathematical background of the inverse problem and optimization scheme. Following, we demonstrate the approach with numerical examples considering different structural geometries subject to various loading conditions. Thereafter, an analysis of the the effects of discretization and imposed constraints is conducted. Finally, discussion and conclusions are provided.

Inverse problem
In this section we propose a regime for optimizing the elastic properties of structural elements using methods rooted in solving inverse problems. We begin by describing inverse problems and how they contrast with classic topology optimization problems. The aim of an inverse problem (for one state) is to find the best parameters ϑ resulting in where F is a linear or non-linear operator and d is a data set. In the case of a discretized elasticity inverse problem (i) F is a numerical forward model using, for example the finite element method or finite difference method and (ii) d is a force vector or, most often, displacement field [71]. To solve such a problem, we aim to minimize a functional ψ, with the basic form where ||·|| represents the appropriate norm for the problem (commonly the L 2 or L 1 norm). It is important to recognize that the structure of the inverse problem is quite different than the structure of a topology optimization functional, which -in a classical sense -aims to minimize the compliance where U and K are the displacement field and stiffness matrix for the geometry of interest, respectively. The fundamental objectives of classical topology optimization (TO) and optimization using inverse methods (IM) may now be directly contrasted: • TO: minimize C given boundary conditions, forces, and design constraints.
• IM: minimize ψ given data set d, operator F , constraints on ϑ, and prior knowledge related to the structure of ϑ.
Some advantages of using inverse methods in elasticity optimization problems are their flexibility to incorporate (a) prior knowledge related to the spatial distribution/structure of ϑ (usually incorporated in regularization techniques), (b) stacked models, i.e. simultaneous solutions for multiple states, and (c) constraints on ϑ, which may result from physically-realistic ranges of a given elastic property or manufacturing limitations. The ability to use (a-c) in determining optimized solutions of elastic properties is the primary motivator for proposing the inverse regime herein.

Problem Statement
The inverse problem is described as follows: given an unloaded structural geometry Ω, boundary information ∂Ω, N target displacement fields with each state's displacement target concatenated in ∆ n = {∆ 1 , ∆ 2 , ∆ 3 , . . . , ∆ N } (in other words, ∆ n is a stacked vector of displacement fields conforming to any design constraints for structural state n), and loading condi- From a more classical viewpoint, we observe that Eq. 5 is equivalent to N systems of concatenated structural equations with each individual state n defined by where K n is a stiffness matrix determined numerically. The concatenated classical model may then be written in a more classical description using: In the following subsection, we will formulate the solution approach to the model problems in Eqs. 5 and 7.

Inverse Approach
We now aim to use the equivalent model problems in Eqs. 5 and 7 to formulate a least-squares (LS) minimization regime to solve the inverse problem. For illustration purposes, we begin by considering only one state, where N = 1 and the problem is unconstrained. In such a case, we aim to determine E 1 by minimizing the functional where || · || represents the L 2 norm, the superscript 1 denotes the degree of n (1, in this case), λ is the Tikhonov regularization parameter [66] which will be detailed later in this subsection, the regularization matrix is given by Γ = λI, and I is the identity matrix.
In general, however, we note that physicallyrealistic (non-complex) solutions are realized when E 1 > 0 -or more generally E n > 0. Moreover, when prior information about the material is available, we may further constrain the problem such that the elasticity modulus is greater than some constant E c known a priori. While still considering N = 1, we may rewrite the functional to be minimized as The extension of Eq. 9 for N = 2 is straightforward, and in the case of equal weighting of cost functions for each state n, is done by summation; i.e. for N = 2 we aim to minimize In the general case, where N may take any positive integer, we aim to find the elastic modulus distributions E n 1 that minimizes the following functional where ω n is the user-defined weight for each state. The inclusion of regularization in Eqs. 8 -11 results from the ill-posed nature of the inverse problem, meaning that standard LS approaches may yield non-unique solutions. The form of regularization chosen here (L 2 regularization) is known to result in smooth solutions [67]. The primary advantage of this regularization technique is its simplicity: the smoothness of solutions is controlled by the magnitude of λ. For further information and discussion on more advanced regularization techniques incorporating structural prior information, we refer the reader to [68,69]. We would like to note that (i) the constraints E n > E c are handled using the interior point method and (ii) here, each state n is equally weighted during the minimization of Ψ such that ∑ N n=1 ω n = 1; however, biased weighting of each state n may also be employed in a straight-forward manner.
To compute the simulated displacement fields 2 U n we employ the finite element method using where q is the total number of unknown displacements, j, i denote the column and row locations, and K −1 ji and F i are often referred to as the compliance matrix and force vector, respectively [70,71]. In this work, we use piece-wise linear triangular elements and neglect the effects of out-of-plane displacements (a plate in plane stress) in generating the stiffness matrix K ji . We would like to comment that the structure of the wellknown stiffness (or compliance) matrices within the FEM equations consisting of F = KU or U = K −1 F are positive definite (when properly constructed), which leads to unique solutions to the forward problem [71]. The ill posedness discussed herein refers to the inverse problem, wherein solutions (in this case E n ) to the inverse problem are highly sensitive to modeling errors and outlier data [68].
To solve the optimization problem we use a Gauss-Newton scheme equipped with a line-search function to compute the step size s k in the parameterized (stacked) solution θ k = θ k−1 + s kθ where θ k is the current estimate andθ is the LS update at iteration k.Within the parameterization, we define notation for the stacked vectors as follows: θ k = [E 1,k , E 2,k , . . . , E N,k ] T and θ k−1 = [E 1,k−1 , E 2,k−1 , . . . , E N,k−1 ] T . For example, E 1,k refers to E n=1,k or "E n=1 at the k th iteration." Specifically, we define the stacked parameterizations and the stacked LS updatē . . .
where H E n and g E n are the Hessian matrix and gradient vector, respectively, computed from the cubic polynomial barrier functions (for additional details regarding the implementation of barrier constraints in a LS framework, like those used herein, we refer the reader to ( [72,73]), and J E n = ∂U n (E n ) ∂E n is the Jacobian. In this work, the Jacobian at each state J E n = ∂U n (E n ) ∂E n is computed at each iteration using the perturbation method with central differencing, where the entries are calculated using where the perturbation p is computed as a function of the machine precision ε using p = 3 ε 2 following [74]. As with many optimization problems, some stopping criteria ϕ must be defined; in this work we found ϕ = (Ψ k − Ψ k−5 )/Ψ k−5 ≤ 10 −3 to be satisfactory. In particular, the use of Ψ k−5 was selected to help ensure that the cost function has reached a stable minimum in the case there are oscillations at later iterations. Lastly, the final optimized distribution of the elasticity modulus E Ψ must be computed. To be consistent with the equal weighting employed in the minimization of Ψ in Eq. 11, we simply compute E Ψ using the mean value for each discretized element in the solution vector θ.

Remark:
We would like to mention that the proposed method is also versatile to optimize other structural and material parameters. For example, one may modify Eqs. 11, 12, 13, 14, and 15 in a straight-forward manner to solve the optimization problem with respect to the Poisson ratio ν, shear modulus G, or plate thickness t. In cases where multiple parameters, for example E and t, are optimized simultaneously, one may incorporate a parameterization regime such as the one proposed in [75,76].

Numerical examples
The following numerical examples are designed to test the proposed approach and demonstrate optimized distributions of the elastic modulus E Ψ in structural elements with N = 2. In the first two examples, we examine the rectangular beam-columns. For the first two examples, (i) the design displacement targets ∆ 1 and ∆ 2 are selected as simple functions in the x, y coordinate plane and (ii) the constraint E c was selected such that the displacement fields were small during the initial iterations. In the third example, we study the effect of mesh size the constraint E c on the inverse solutions for an irregular geometry; we also consider the numerical behavior during minimization of the cost function. In the examples, we neglect out of plane deformations, i.e. the structural elements are in plane stress. For simplicity, we use dimensionless units for the deformations, forces, elastic coefficients, and coordinate system.
In solving each problem, we select λ = 10 −3 and make no further attempt to optimize the regularization parameter during iterations. Lastly, all problems were solved using a MATLAB implementation of the algorithm described in section 2 and the choice of the domain meshing was selected based on available computational resources.

Example 1: non-slender beam column
In this example, we consider a cantilevered beam clamped at y = 0 with two design displacement targets, a homogeneous Poisson ratio ν = 0.35, and plate thickness t = 0.2. The first design displacement target takes the form of a cubic polynomial where the mesh points of the unloaded beam are displaced by the following function ∆ 1 (x, y) = [0.005y 3 − 0.07y 2 + 0.2y, −], where the y coordinate refers to the vertical coordinates along the beam (note that the y-coordinates were intentionally not translated). The second target is given by a pure axial translation of the form ∆ 2 (x, y) = [−, 0.01y]. The design displacement targets and loading for each criterion are shown in Fig. 1 plotted atop the displaced triangulation consisting of 800 elements. Moreover, in determining E Ψ , we utilized the constraint E c = 30. The solution for this example is shown in Fig.  2, with optimized E Ψ plotted atop displaced configurations generated from simulated displacement fields F 4 ), which correspond to displacement targets ∆ 1 and ∆ 2 , respectively. Results shown in Fig. 2 show a highly inhomogenous distribution of the elasticity modulus. This is a satisfactory and expected result, given (a) the non-physical nature of ∆ 1 , where the y-components of the field were not translated as would be observed in physically-realistic bending case (i.e. ∆ 1 did not consider rotations) and (b) the applied forces are not congruent with the displacement fields prescribed in ∆ 1 and ∆ 2 . It is therefore apparent why areas of the beam have been stiffened and softened to compensate for discrepancies between the forces and prescribed displacement fields. It is interesting to note that near the location of the largest applied point force, F 2 , there is an obvious local increase in E Ψ . This results from local data mismatch between the smooth field ∆ 1 and the locally-perturbed field U 1 . The local increase in stiffness near F 2 is therefore consistent with this realization. Moreover, upon comparing Figs. 1 and 2 we observe that U 1 and ∆ 1 are not a perfect match -despite meeting the minimization criterion ϕ. On the other hand, we observe that U 2 and ∆ 2 compare more favorably. Central reasons for these observations are (i) the previously noted unrealistic nature of ∆ 1 (in contrast to ∆ 2 , which is a realistic description of an axial displacement field), (ii) the ill-posed nature of the inverse problem, (iii) the increased non-uniqueness in stacked parameterizations of inverse problems (also noted in [75,76], (iv) the use the L 2 regularization, an uninformative prior model, (v) the regularization parameter λ was not optimized at each iteration, and (vi) the effects of weighting the solution E Ψ .

Example 2: slender beam column
In this example, we consider a slender beam with t = 0.1 and ν = 0.35. The beam is clamped at y = 0 and pinned at the coordinates (x, y) = [1,10] without the presence of axial loads. The first displacement target is quadratic and has the form ∆ 1 (x, y) = [−0.015y 2 + 0.15y, −]. The second target is a cubic function given by ∆ 1 (x, y) = [0.005y 3 − 0.08y 2 + 0.25y, −]. No translation in the y-direction is prescribed. The target displacement fields and loading for each criterion are shown in Fig. 3 plotted atop the displaced triangulation consisting of 980 elements. The constraint was chosen as E c = 100. The solution for this example is shown in Fig.  4, with optimized E Ψ plotted atop displaced configurations generated from simulated displacement fields F 2 , F 3 ), which correspond to displacement targets ∆ 1 and ∆ 2 , respectively. We again observe a highly inhomogenous distribution of the solution E Ψ , especially in the areas of high curvature. As in Example 1, the local stiffening and softening observed in E Ψ is largely affected by the local discrepancies between the simulated displacement fields and target displacement fields. Of note, significant stiffening is observed near the pinned connection at (x, y) = [1, 10] and the negative moment region at the clamped base as the inverse model is compensating for localized effects resulting from high deformations gradients. Relative to the beam in Example 1, the beam in this example is markedly more slender. This led to a significant differences in the distribution of the elastic modulus in the solution. Upon observation of Fig. 4, it is apparent that the variability of E Ψ is much higher than the variability observed in Fig 2. This is most plausibly due to the lower shear stresses occurring in the more slender beam. In the case of slender beams, differences in the simulated and target displacement fields are most significantly compensated by changes in the bending stiffness near the top and bottom of the beam cross sections. On the other hand, solutions for deeper beams are more dependent on effects near the neutral axis of the member due to shear effects. Bearing this realization in mind, is important to note that, while differences in the aspect ratio of the member do have important implications on the global distribution of E Ψ , the effects of external loads and boundary conditions also have significant local implications on E Ψ .
3.3 Example 3: Irregular member considering the effects of mesh size and constraints In the final example, we consider an irregular member with t = 0.25 and ν = 0.35. The member is clamped at y = 0. The first target displacement is of the form ∆ 1 = [0.005y 2 , −0.005x 2 ] and the second target displacement is a linear translation of the form ∆ 2 = [0.02y, −]. Two discretizations are considered, the first mesh, denoted the coarse mesh, consists of 67 triangular elements. The second discretization considered has 247 elements, and is denoted the refined mesh. The target displacement fields and loading for each criterion are shown in Fig. 5 plotted atop the displaced triangulations.
The solutions for this example are shown in Fig.  6, with optimized E Ψ plotted atop displaced configurations generated from simulated displacement fields , which correspond to displacement targets ∆ 1 and ∆ 2 , respectively. We observe that, in general, the obtained solutions have similar distributions of E Ψ . Indeed, we clearly observe, in all images, that the column portion of the member is significantly more stiffened relative to the horizontal member due to the bending moments resulting from F 2 . On the other hand, since the horizontal member is in a state of axial tension in the second state, the optimized stiffness for the horizontal component is considerably lower than the vertical component. This is therefore an expected result.
In comparing the results obtained using different discretizations, we note two primary differences in the solutions: the distribution of the elastic modulus at the base of the column and near the location of F 1 . In solutions obtained using the coarse mesh, there is significant stiffening near the base and nodal force. This observation results from two realizations, (a) that coarser meshes are known to decrease model accuracy during stiffness integration, in many cases increasing stiffness [70] and (b) coarser meshes decrease the degrees of freedom for solutions (E Ψ ) in the inverse problem. Based on (a), it is no surprise that the increased element and global stiffness resulting from coarse forward modeling is reflected in the inverse solutions. With regards to point (b), it is apparent that by reducing the degrees of freedom in E Ψ , the coarse solutions are smoother since small localized fluctuations of E Ψ are not possible. Solutions obtained with the fine mesh, on the other hand, do have small fluctuations which may be viewed as artifacts resulting from the non-uniqueness of the inverse problem and the non-optimized regularization parameter [68].
For results reported in this section, the effect of the lower constraint E c on the optimized solutions of E Ψ is more subtle than the effect of discretization. The effect of E c was similar for solutions obtained using coarse and fine meshes: the results were smoother for E c > 50 than for E c > 0. The reason for this is that the available solutions 3 of E Ψ are more restricted and closer to the average optimal values of E Ψ .
Of particular interest in optimization problems is the minimization rate of the cost function, which sheds light on the computational efficiency of the regime. Here, we compare the minimization rate for both discretizations and consider the effect E c has on the final solution once the stopping criteria is reached. The minimization curves for the cost function Ψ as a function of the iteration number k are shown in Fig.7. The minimization curves shown in Fig 7 support the earlier statements regarding the behavior of the optimization problem with respect to the effects of the constraint E c and the mesh refinement; this will be discussed in the following. While the minimization rates of Ψ for both the discretizations were similar, the refined meshes reached a lower minimum for both E c > 0 and E c > 50. Moreover, for both the coarse and fine meshes, E c > 0 resulted in lower minimum values of Ψ. These results indicate that by increasing the degrees of freedom available for the inverse solution (via decreasing E c and increasing the spatial degrees of freedom), the optimized solution for E Ψ reached a more globally optimized solution.

Discussion
In this section, we examined optimized solutions of the elastic modulus E Ψ for 2D rectangular beams and an irregularly shaped beam-column considering different loading, boundary conditions, target displacement fields, and aspect ratios using the proposed regime. In the analysis of the irregularly-shaped beam column, the effects of the lower constraint E c and discretization were investigated. As a whole, the proposed opti-mization regime and the solutions obtained were determined to be physically and intuitively consistent. We noted multiple factors contributing the inhomgeneity, variability, and distribution of the optimized solutions. Some of the important contributors in solution of E Ψ included (in no specific order): 1. Global effects related to the geometry of structural elements. 2. Local effects resulting from external forces and imposed boundary conditions. 3. The (possibly) non-physical nature of prescribed displacement fields ∆ n . 4. Effects from not optimizing the regularization parameter λ. 5. Effects from weighting E Ψ . 6. Local and global effects due to mesh size. 7. Global effects related to the lower constraint E c .
It is tempting to draw general conclusions as to which variable contributes most to the optimized solution E Ψ or to which parameter the inverse problem is most sensitive. However, when we consider the coupled non-uniqueness of the inverse problem, number of structural states stacked, complexity of boundary/loading conditions, and prior information; we may surmise that solutions of E Ψ are highly application specific.
We would like to remark that this article has focused primarily on the numerical regimes, minimization techniques, and interpretation of the results in idealized cases. However, it should be noted that current manufacturing techniques used to generate spatial inhomogeneity in elastic properties such as the Young's modulus or Poisson's ratio have limitations. In such cases, we may expect manufacturing processes to have band-limited ranges and constraints. These constraints may, indeed, be used within the constrained optimization scheme by incorporating upper and lower barrier functions at the extremes of the manufacturing bounds. As results of this article have shown, the use of constraints decreases the degrees of freedom and may result in a less optimized solution relative to a problem contrained to be only physically realistic (e.g. E > 0).

Conclusions
In this article, we proposed a constrained leastsquares approach for computationally optimizing the distribution of the elastic modulus for N states. Each state was defined by different loading conditions and target displacement fields that may be prescribed by the designer. A numerical study was conducted to test the approach, where 2D structural elements with different aspect ratios, geometries, loading/boundary conditions, and target displacement fields were analyzed. The approach was tested for N = 2, where results were demonstrated to be physically-realistic and intuitively consistent. It was noted that optimized distributions of the elastic modulus were affected by important factors, such as the structural geometry, prior information, solution weighting, physically-realistic nature of the target displacement fields, imposed constraints, discretization, and overall non-uniqueness of the inverse problem.