Optimization of a lattice spring model with elastoplastic conducting springs: A case study

To model multi-functional materials with plasticity, this paper proposes to employ the so-called lattice spring model where the sweeping process theory [15; 6] recently provided algebraic formulas [5] for the response force and where the springs can simultaneously be viewed as conductors. Efficient methods that map actual materials to lattice spring models are developed in e.g. Chen et al [3] and Kale and Ostoja-Starzewski [10].

The analysis in this paper concerns the lattice spring model of Fig. 1, because this is the minimal irreducible model that requires a 3D (i.e. non-planar) sweeping process framework (Rachinskiy [18]). Therefore, we expect that all interesting non-trivial optimization phenomena are to be captured by the model of Fig. 1. In other words, we expect that model of Fig. 1 is an efficient benchmark to assess the capabilities of the sweeping process framework in the optimization of most typical multi-functional properties of lattice-spring models of general form. Indeed, using the model of Fig. 1 we conjectured a number of qualitative properties of universe nature (see Conclusions section). On the other hand, we included quotes from Wolfram Mathematica code and a computational guide (Appendix 8) that allow to generalize simulations of this paper to arbitrary lattice-spring models with 1-dimensional nodes (as long as the resistance of the network is computable via the known methods of e.g. [4; 9]). Further extension of the sweeping process framework to the case of nodes of arbitrary dimension and hyperuniform network materials is available in [8] whose formulas can also be used for optimization of multi-functional properties along the strategy of the current paper.

In what follows, $c_{i}$ denotes the elastic limit of spring $c_{i}$, so that the relaxed length of the spring changes when the stress attempts to exceed the interval $[-c_{i},c_{i}].$ For lattice spring models coming from fiber-reinforced composites, the sum $c_{1}+c_{2}+c_{3}+c_{4}+c_{5}$ can be viewed as the fabrication cost of the material because the cost of such composites is mainly determined by fabrication of the nano-carbon fibers [14]. When the displacement-controlled loading (applied to nodes ① and ④) stretches the system gradually, the response force of the system will increase up to a certain terminal value $F$ that quantifies the strength of the material. Viewing Fig. 1 as an electric circuit, and using $R$ and $G$ to denote the resistance and conductance between nodes ① and ④, the present paper addresses optimization of the functionals

Using Wolfram Mathematica we make a variety of qualitative observations (listed in Conclusions section) that could motivate a good deal of new theoretical research in optimization theory. Some of these observations are uniqueness of optimal values of parameters $(c_{1},...,c_{5})$ in the optimization of $F_{R}$ and $F_{G}$, non-uniqueness of $(c_{1},...,c_{5})$ in the optimization of $C$, existence of linear thresholds in the space of parameters that separate topologically different patterns of optimal $(c_{1},...,c_{5}).$

The paper is organized as follows. In the next section of the paper we introduce the ingredients $F$, $R$, and $G$ used in the definition of the functionals $F_{R}$ and $F_{G}$. In particular, we introduce two domains (6) and (8) for the parameters $c_{1},...,c_{5}$ that this paper restricts the analysis to (taken from Gudoshnikov et al [5]). Section 3 discusses the numerical methods used for optimization of the functionals $F_{R}$, $F_{G}$, $C$, along with computational complexity of the method of our choice. Section 4 is devoted to maximization of the functionals $F_{R}$ and $F_{G}$ with a constraint on the cost $C.$ In this section we first determine the maximal value $F_{max}$ of $F$ (Section 4.1) that can be reached by the system of Fig. 1 with the chosen constraint on the cost and then maximize $F_{R}$ and $F_{G}$ (Sections 4.2 and 4.3 respectively) restricting on the domain of $c_{1},...,c_{5}$ for which the value of $F$ doesn’t go below $F_{max}$ by more than $25\%$ (constraint (13)). In other words, we forbid improved multi-functional properties of a lattice-spring system to compromise the strength of the system by more than $25\%.$ Minimization of the cost $C$ under constraint on $F_{R}$ and $F_{G}$ (from below) is addressed in Sections 5.1 and 5.2. Sections 5.1 and 5.2 are restricted on the domain of $c_{1},...,c_{5}$ for which the value of $F$ doesn’t go below $F_{max}$ too. The conclusions from numerical simulations of Sections 4 and 5 are formulated as Propositions 4.2, 4.3, 5.1, 5.2. Additionally, Section 4.1 quotes the Mathematica’s command we use, which clarifies that all numerical computations can be executed analytically by hand over the methods of linear programming. The conclusion section interprets the numerical results of Sections 4 and 5.1 from more general theoretic perspectives and draws a list of analytic questions inspired by Propositions 4.2, 4.3, 5.1, 5.2. Two appendixes that compute the resistance of the circuit of Fig. 1 (Appendix 7) and outline the sweeping process framework of [5] based conclusions about plastic regimes in elastoplastic system of Fig. 1 (Appendix 8) conclude the paper. Wolfram Mathematica Notebook is attached as supplementary material.

Resistance of the System of Springs. Assuming that the elastic limits of the 5 springs of Fig. 1 are given by $[-c_{i},c_{i}]$, $i\in\overline{1,5}$, the resistance $R_{i}$ of spring $i$ will be taken as

The relation (2) is inspired by typical design mechanisms used in composite materials. For example, one possible realization of the spring can be achieved by fiber-reinforced polymer matrix composites [14]. In this system, both the mechanical strength (i.e. $c_{i}$) and conductivity (i.e. $1/R_{i}$) are determined by the fiber phase and thus, are positively correlated [21]. With the values of resistances of springs given by $R_{i}$, the resultant resistance between nodes ① and ④ computes as (see Appendix 7)

and conductance is computed as

Terminal Response Force. According to the sweeping process framework by Gudoshnikov et al [5] (see Appendix 8 for details), gradual stretching of the lattice of Fig. 1 via the displacement-controlled loading that is applied between nodes ① and ④ increases the response force of the lattice up to a certain terminal value whose formula takes different forms depending on the values of $c_{1},...,c_{5}$. These formulas can be summarized as follows

where we restrict ourselves to parameters that ensure the greatest number of springs reaches the plastic mode. The later property corresponds to the most uniform distribution of plastic deformation across the lattice that reduces the risk of crack initialization. Specifically, if inequalities (6) hold in a strict sense then springs 1, 3, 5 reach plastic deformation as the magnitude of the displacement-controlled loading crosses a threshold (i.e. springs 1, 3, 5 get to the plastic mode terminally, see [5] for explicit formulas), while if condition (8) holds in a strict sense then those are springs 2, 3, 4 that get to plastic mode terminally¹¹1Note, the non-strict inequalities in (8.4) and (8.5) of [5] is a typo. There should be strict inequalities because (8.4) and (8.5) come from (7.10) whose inequalities are strict.. The inequalities (6) and (8) will be referred to as feasibility conditions. In what follows, we use feasibility conditions as optimization constraints, so that we will usually get optimal $c_{i}$ on the boundary of the domains (6) and (8) and strict inequalities in (6) and (8) won’t hold for our optimal $c_{i}$. Arbitrary small perturbation will be needed to move such optimal $c_{i}$ to the interior of the domain (6) or (8) to make the plastic behavior of the model predictable according to [5].

We analyze the results for the optimization problem based on two different optimization methods for constrained global optimization available within Wolfram Mathematica. We use NMaximize and NMinimize which are global optimization algorithms for solving problems numerically. Within those we have other optimization methods, the ones that we consider for our problem are ”RandomSearch” and ”Differential Evolution”. Differential evolution is a simple stochastic function minimizer. It starts [23; 24] with a population of candidate solutions randomly distributed in the search space. At each iteration, new candidate solutions are generated by combining and mutating existing solutions according to certain rules. These new solutions are then evaluated, and the best ones are retained for the next generation. Random search on the other hand works [23] by generating a population of random starting points and uses a local optimization method from each of the starting points to converge to a local minimum. When the results from the two methods do not coincide, Differential Evolution was the one showing the correct extremum except of a negligible number of exceptions that we furter verified using Phyton SLSQP algorithm. We therefore conclude that Differential Evolution is a reliable method for the optimization problem considered in this paper.

According to [17], the computational complexity of the method of Differential Evolution for large networks is determined by the dimensions of the search space, the population size, and the number of iterations. Specifically, the general computational complexity for Differential Evolution is given by [17]

under assumption that the computational complexity of evaluating the objective function grows at most linearly with the dimension of the search space (which is the case for the objective functions of the present paper because $C$ is linear by construction, $F$ is linear by [6, formula (3.23)], and computational time for $R$ scales close to linear if using the algorithm [13], see [13, table I]). Formula (10) says that the number of operations required by Differential Evolution is proportional to the product of the problem’s dimensionality ($n$), the population size ($N_{p}$), and the number of iterations ($G_{max}$). The population size $N_{p}$ is usually a fixed constant between $5n$ and $10n$ [20]. The number of iterations $G_{max}$ depends on the desired accuracy and can be viewed as another fixed constant independent of $n$ and $N_{p}$. Therefore, if the population size is proportional to $n$, then computational complexity of the method of Differential Evolution scales quadratically with respect to the dimension of the problem.

In this section we search for the parameters $c_{i}$ that optimize $M_{R}$ and $M_{G}$ under assumption that

Since the optimal values of the functionals $F_{R}$ and $F_{G}$ vary in the intervals [0.1,1.2] and [0.1,1.5] when the cost of the material is constrained (Sections 4.2 and 4.3), we will pick a constant inside both of the intervals as a constraint (lower bound) for the functionals $F_{R}$ and $F_{G}$ when optimizing $C$. Specifically, we will consider

and

in Sections 5.1 and 5.2 to follow.

The conclusions obtained with Wolfram Mathematica NMaximize command are summarized in Propositions 4.2, 4.3, 5.1, 5.2 and motivate several tasks towards rigorous mathematical proofs that we group below according to the above stated Problems A, B, C, and D. We recall that we refer to an optimal value $(c_{1},...,c_{5})$ as unique, if $(c_{1},...,c_{5})$ is unique in one of the domains (6) and (8), see Remark 1.

A1: Prove uniqueness of the optimal value $(c_{1},...,c_{5})$ for each value of $(k_{1},k_{2}).$

A2: Prove that the space of parameters $(k_{1},k_{2})$ splits into 3 groups so that the optimal values $(c_{1},...,c_{5})$ in each group coincide.

A3: Examine whether or not one of the 3 groups from A1 always forms a straight line that separates the other 2 groups. At the very least, prove that the above-mentioned 2 groups are always separated by a straight line $L$. Find the equation of $L$, which is expected to be independent of the chosen domain of parameters $(c_{1},...,c_{5})$ (domains (6) and (8) in this paper).

A4: Prove that the optimal values $(c_{1},...,c_{5})$ below $L$ belong to the boundaries of both domains (6) and (8) simultaneously. Prove that the optimal values $(c_{1},...,c_{5})$ above $L$ and on $L$ do not belong to the boundaries of (6) and (8) simultaneously.

B1: Prove that optimal values $(c_{1},...,c_{5})$ are unique, independent of parameters $(k_{1},k_{2}),$ and belong to the boundaries of both domains (6) and (8) simultaneously.

C1: Prove the existence of a linear threshold $L_{2}$ in the space of parameters $(k_{1},k_{2})$ on so that the optimal values $(c_{1},...,c_{5})$ above $L_{2}$ coincide. Find the equation of $L_{2}$, which is expected to depend on the chosen domain of parameters $(c_{1},...,c_{5})$.

C2: For the values of $(k_{1},k_{2})$ below $L_{2}$, prove the existence of lines (presumably parallel to $L_{2}$) so that the optimal values of $(c_{1},...,c_{5})$ are constant alone these lines.

C3: Prove the existence of another linear threshold $L_{1}$ in the space of parameters $(k_{1},k_{2})$ (below $L_{2}$) such that $C$ keeps a constant value above $L_{1}$ and $C$ is larger than this constant value below $L_{1}.$ Find the equation of $L_{1}$, which (in contrast with $L_{2}$) is expected to be independent on the chosen domain of parameters $(c_{1},...,c_{5}$ and expected to be parallel to $L_{2}$.

C4: Prove that crossing the threshold $L_{1}$ corresponds to $F_{R}$ reaching its chosen bound (in this paper, crossing $L_{1}$ from above to below, $F_{R}$ reaches the lower bound of $0.5$).

C5: Prove that $F$ equals its chosen bound (0.75 in this paper) for all $(k_{1},k_{2})$ except for isolated special values. Determine the nature of these special values.

C6: Prove that, across all values of $(k_{1},k_{2})\in[k_{1,min},k_{1,max}]\times[k_{2,min},k_{2,max}]$, the maximum of the minimal values of $C$ is attained at some $(k_{1,*},k_{2,*})\not=(k_{1,min},k_{2,min}).$ Moreover, such a value of $(k_{1,*},k_{2,*})$ is unique.

C7: Prove the existence of a linear threshold such that the optimal values $(c_{1},...,c_{5})$ are unique above the threshold and non-unique below the threshold.

C8: Prove that the optimal values $(c_{1},...,c_{5})$ above $L_{1}$ belong to the boundaries of both domains (6) and (8) simultaneously. Prove that the optimal values $(c_{1},...,c_{5})$ below $L_{1}$ do not belong to the boundaries of (6) and (8) simultaneously.

D1: Prove the existence of a linear threshold $L$ in the space of parameters $(k_{1},k_{2})$ such that $C$ keeps a constant value above $L$ and $C$ is larger than this constant value below $L.$ Find the equation of $L$, which is expected to be independent of the chosen domain of parameters $(c_{1},...,c_{5})$.

D2: Prove that crossing the threshold $L$ corresponds to $F$ switching from $F>const$ to $F=const$ ($const=0.75$ in this paper) and switching from $F_{G}=const$ to $F_{G}>const$ ($const=0.5$), if moving $(k_{1},k_{2})$ from under $L$ to above $L$.

D3: For the values of $(k_{1},k_{2})$ below $L$, prove the existence of lines along where the optimal values of $(c_{1},...,c_{5})$ are constant.

D4: Prove uniqueness of the optimal values $(c_{1},...,c_{5})$ for every $(k_{1},k_{2})$ potentially disproving the existence of those special pairs $(k_{1},k_{2})$ (such as $(k_{1},k_{2})=(0.2,0.6)$ and $(k_{1},k_{2})=(0.2,0.4)$ in Fig. 11) where nonuniqueness is detected by Wolfram Mathematica (or explain the nature of $(k_{1},k_{2})$ where nonuniquness of the optimal values takes places).

Simulations of this paper are executed under assumption that the plastic mode is terminally reached by a maximally possible number of springs. Dropping this assumption might be of interest in applications outside of materials science. This would require extension of our numerical analysis from domains (6) and (8) to all positive values of $c_{1},...,c_{5}$, that might stimulate further development of the theory of [5] and further derivation of the corresponding terminal forces (6) and (8) (that are not yet available for all positive values of $c_{1},...,c_{5}$).

We hope that simulation results reported and theoretical questions motivated by these simulation results will attract interest of future researchers.