## 2 Theory  
### Optimization Formulae  
The objective function in the optimization problem is constructed using the relationship between the prescribed values on the F-D plane and the estimated values from the layouts in the optimization process. The displacements at the defined forces have to be calculated for each iteration of the optimization process. Structures show nonlinear behavior when they undergo large deformation and usually have inherent monotonous softening or hardening nonlinearities with respect to external forces. For instance, a thin rectangular plate exhibits monotonous hardening nonlinearity. In this work, we focus on monotonic softening and hardening nonlinearities, which allows us to use the Newton-Raphson method to calculate the displacement field. The solution does not converge when using a small number of target points on the F-D plane and a rectangular thin plate as the initial design model. This occurs because of the difference between the inherent nonlinearity in the initial design model and the desired nonlinearity. Further descriptions can be found in Section 2.3. To help with the nonlinearity adaptation during the optimization process, the definable slopes at the target points are used as constraints in the optimization problem. Since we are dealing with a density-based topology optimization problem, an individual design variable \(\rho\) can vary continuously in the interval \([0,1]\), where \(0\) and \(1\) imply the absence and presence of material, respectively. On the basis of these assumptions, the optimization problem can be defined as follows:  
\[ \begin{array}{ll}\min\limits_{0<\rho<1}&:c(\rho)=\frac{N_{t}}{i}w_{i}^{c} \Big{[}\mathbf{I}^{T}\mathbf{u}(f_{i},\rho)-u_{i}^{*}\Big{]}^{2}\\ s.t.&:\ \mathbf{r}\ =\ 0\\ &:\ \theta_{i}^{\min}\leq\theta_{i}^{c}(f_{i})\leq\theta_{i}^{\max}\,\ \ i\ =\ 1,\ N_{t}\\ &:\ \frac{(\overline{\rho})}{\overline{V_{0}}}\leq V_{\max}\end{array}  \tag{1}\]  
\(c\) is the objective function, and \(\rho\) is a vector of design variables, \(N_{t}\) is the number of defined target points on the F-D curve, \(f_{i}\) is the amplitude of the applied \(i\)th external force, \(\mathbf{u}\) is the displacement vector in global coordinate, \(u_{i}^{*}\) is the prescribed value of the displacement with respect to the \(i\)th external force, and \(w_{i}^{c}\) is the weighting function. \(\mathbf{I}\) is the position vector, in which a unit value is assigned at the corresponding degree of freedom where displacements are identified, whereas the other components are set to zero. The superscript T indicates the transpose of a vector. \(\theta_{i}^{c}\) is the slope that can be defined on the F-D plane for the \(i\)th external force. Details about the slope and its derivative are given in Section 2.3. \(\theta_{i}^{\min}\) and \(\theta_{i}^{\max}\) are the upper and lower bounds of the slope constraint, \(\overline{\rho}\) is a vector of physical density, and \(V_{0}\) is the volume when the entire design domain is occupied by solid elements. \(V_{\max}\) is the upper bound of the volume fraction constraint, and \(\mathbf{r}\) is the residual force vector of the structure equilibrium. The residual force vector is iteratively estimated in the nonlinear analysis. Large deformation with a small strain state is assumed in this study. This means that there is a linear relationship between the Green-Lagrangian strain tensor and the second Piola-Kirschhoff stress tensor \(\mathbf{S}\):  
\[\mathbf{S}=\mathbf{D}:\mathbf{E} \tag{2}\]  
The fourth-order tensor \(\mathbf{D}\) is the elastic tensor, and \(\mathbf{E}\) is the Green-Lagrangian strain tensor, which is represented as \(\mathbf{E}=\Big{(}\mathbf{F}^{T}\mathbf{F}\ -\ \mathbf{I}\Big{)}/2\). \(\mathbf{I}\) is the identify matrix. \(\mathbf{F}\) is the deformation gradient tensor defined as \(\mathbf{F}\ =\ \mathbf{I}\ +\ \partial\mathbf{u}/\partial\mathbf{X}\). \(\mathbf{X}\) denotes the material coordinate in the undeformed geometry. The total Lagrangian finite element formulation is used.