Sage Journals: Discover world-class research

Abstract

Aiming at the problem that many algorithms could not effectively balance the global search ability and local search ability, a new optimization algorithm is proposed. Inspired by Finite Element Analysis (FEA) approach, a relationship of mapping between Finite Element Analysis approach and a population-based optimization algorithm is constructed through comparing the similarities and differences of FEA node and ideal particle. In algorithm framework, the stiffness coefficient corresponds to a user-defined function of the value of an objective function to be optimized, and the node forces among individuals are defined and an attraction-repulsion rule is established among them. The FEA approach that can simulate multi- states of matter is adopted to balance the global search ability and local search ability in the novel optimization algorithm. A theoretical analysis is made for algorithm parallelism. The conditions for convergence are deduced through analyzing the algorithm based on discrete-time linear system theory. In addition, the performance of the algorithm is compared with PSO for five states which include free state, diffusion state, solid state, entirely solid state, synthesis state. The simulation results of six benchmark functions show that the algorithm is effective. The algorithm supplies a new method to solve optimization problem.

Keywords

Finite element analysis optimization multi-states of matter parallelism convergence

Introduction

Inspired by the laws of nature and swarm intelligence in biology society, heuristic algorithm applies different strategies to search the solution space. Genetic algorithms mimic Darwinian natural selection,¹ where “fitness” selects individuals for inheritance, selection, crossover, and adaptive mutation. Particle Swarm Optimization (PSO) which is based on the swarming behavior of birds seeking food,² is the earliest example of a Swarm Intelligence algorithm. Particle Swarm Optimization (PSO) algorithm has existed premature convergence for multimodal search problems, and PSO has undergone numerous improvements.^3–4 Some nature-inspired heuristics analogize the laws of physics, as do CFO, EM and APO. CFO searches for extrema by “flying” through the decision space a group of “probes” whose trajectories are deterministically governed by the gravitational analogy.^5–7 EM is inspired Coulomb’s Force Law for electrostatic charges.⁸ A framework of artificial physics optimization (APO) algorithm is constructed, in which a swarm of particles are treated as physical individuals seeking for a global optimum in the problem space driven by virtual force.⁹

These heuristics have their applying ranges, advantages and disadvantages respectively. Some heuristics are easily trapped in the local optima of highly multimodal or “deceptive” functions. Therefore, a lot of problems about heuristic algorithm remain to study further. Finite Element Analysis¹⁰ is one of the important foundational theories of digital technology. FEA has the features of simulating multi-states of mater, which has reference value for innovation of intelligence optimization algorithm. The primary principle of FEA is that a complex system can be divided into many comparative small and simple elements which facilitate understanding complex systems, so that a thoughtful understanding is obtained by combining the elements.¹¹ In many fields finite element method is used for simulating multi-states of matter,^12–13 including solid mechanics, flow simulation, and gas flow field.

Aiming at the problem that many algorithms using a single motion rule could not effectively balance the global search ability and local search ability, a novel optimization algorithm based on Finite Element Analysis which is easy to simulate multi-states of matter and to implement parallel computing is proposed. Compared with the PSO and EM, the proposed algorithm is inspired by different laws of physics. The novel algorithm has advantages of effective motion which is affected by all the other particle, while PSO only flies to the historical best particle and the current best particle. The novel algorithm which has multi- states of matter is better at fine search than EM. Therefore, compared with the PSO and EM, the proposed algorithm has advantages of higher diversity and search efficiency. The paper constructs algorithm framework and verify the effectiveness of the algorithm through experiment. A theoretical analysis is made for algorithm parallelism at main steps, and the implementation of parallelism is a further study work.

The remainder of this paper is organized as follows: Section 2 introduces the FEA approach. Section 3 describes a relationship of mapping between FEA and a population-based optimization algorithm. Section 4 describes FEAO framework and analyzes theoretical parallelism of FEAO algorithm. Section 5 analyzes FEAO convergence. Section 6 compares the performance of FEAO five states and PSO against six recognized benchmark functions. Section 7 presents conclusions and suggestions for future work.

Finite element analysis

The typical finite element analysis goes through roughly four steps.¹⁴ With one-dimension bar structure taken as an object for study, an example is given to intuitively illustrate the process of FEA in Figure 1. The structure consists of n bar components of different geometric sizes, and the junctions of the structure is under concentrated forces. The relevant parameters of material characteristics and structure size: $E_{i}$ (i = 1,2…n) —— the elastic modulus, $A_{i}$ (i = 1,2…n) —— the cross sectional area of the bar, $\frac{E_{i} A_{i}}{l_{i}}$ —— stiffness coefficient of element, $u_{i}$ (i = 1,2…n + 1) —— the displacement of node i, $P_{i}^{(i)} — — the forces on node i .$

(1) Dividing elements and numbering node

Figure 1.

Nodes and elements number of 1D bar structure.

At the junction of structure, n + 1 nodes and n elements are divided. The element and node numbers are illustrated in Figure 1. Each element is separated, and the displacement and external force of each node is indicated.

(2) Element stiffness equation

Stiffness equation (1) of element i (i = 1,2…n) is described as follows:

[\begin{matrix} \frac{E_{i} A_{i}}{l_{i}} & - \frac{E_{i} A_{i}}{l_{i}} \\ - \frac{E_{i} A_{i}}{l_{i}} & \frac{E_{i} A_{i}}{l_{i}} \end{matrix}] [\begin{matrix} u_{i} \\ u_{i + 1} \end{matrix}] = [\begin{matrix} P_{i}^{(i)} \\ P_{i + 1}^{(i)} \end{matrix}]

(1)

(3) Assembling the stiffness matrix and establishing global stiffness equation (2).

\begin{matrix} [\begin{matrix} u_{1} & u_{2} & u_{3} \dots & u_{4} \\ ↓ & ↓ & ↓ & ↓ \\ \frac{E_{1} A_{1}}{l_{1}} & - \frac{E_{1} A_{1}}{l_{1}} & 0 & 0 \\ - \frac{E_{1} A_{1}}{l_{1}} & \frac{E_{1} A_{1}}{l_{1}} + \frac{E_{2} A_{2}}{l_{2}} & - \frac{E_{2} A_{2}}{l_{2}} & 0 \\ 0 & - \frac{E_{2} A_{2}}{l_{2}} & \begin{matrix} \frac{E_{2} A_{2}}{l_{2}} + \frac{E_{3} A_{3}}{l_{3}} \dots \\ \dots ⋱ \end{matrix} & - \frac{E_{n} A_{n}}{l_{n}} \\ 0 & 0 & - \frac{E_{n} A_{n}}{l_{n}} & \frac{E_{n} A_{n}}{l_{3}} \end{matrix}] [\begin{matrix} u_{1} \\ u_{2} \\ u_{3} \\ ⋮ \\ u_{n + 1} \end{matrix}] \\ = [\begin{matrix} P_{1}^{(1)} \\ P_{2}^{(1)} + P_{2}^{(2)} \\ P_{3}^{(2)} + P_{3}^{(3)} \\ ⋮ \\ P_{n + 1}^{(n)} \end{matrix}] \end{matrix}

(2)

(4) Handling boundary conditions and solving global stiffness equation

With the boundary condition $u_{n + 1}$ = 0, (n + 1)th row and (n + 1)th column of global stiffness equation are both eliminated. Global stiffness equation is solved, and subsequently node displacement ${[u_{1} u_{2} u_{3 \dots} 0]}^{T}$ is obtained.

A relationship of mapping between FEA and a population-based optimization algorithm

To address the optimization problem, the solutions sampled from the feasible region of the problem are treated as particle which have properties, such as fitness, velocity, position. With intelligent search strategy, the optimal solution is found by the way of successive iteration velocity and position. The nodes which forces are exerted on will deform appropriately depending on stiffness equations which consist of the forces and stiffness matrix. While nodes of FEA corresponds to particle of optimization algorithm, the paper constructs the relationship between node’s stiffness and the fitness of an objective function to be optimized. In addition, force rule of FEA corresponds to search strategy of optimization algorithm. As above, a relationship of mapping between Finite Element Analysis approach and a population-based optimization algorithm is constructed, as showed in Figure 2.

Figure 2.

The mapping of FEA to optimization algorithm.

Matter can exist as a solid, liquid, or gas. Finite element method is used for simulating multi-states of matter. Adopting FEA simulation capacity, three node states which correspond to three different motion rules are designed for the novel algorithm. By setting the property of element stiffness matrix the nodes enter one of three states which include solid state, diffusion state, and free state. Three main states of FEAO algorithm draw analogies from three motion rules which effectively balance the global search ability and local search ability in optimization algorithm. (1) Solid state: solid state draws analogies from this method which can disperse the particles and improve the diversity of the population in order to avoid premature convergence. According to solid mechanics, atoms of solid states are not easily to move. Strategy of solid state turn the excessive-gathering particles among population into solid state which avoid excessive-gathering particles to get too close. The implement method of solid state is that particular negative values are assigned to off- diagonal elements of stiffness matrix. (2) Diffusion state: Diffusion state is the exact opposite of solid state. In diffusion state the node motion is motivated. The implement method of diffusion state is that particular positive values are assigned to off-diagonal elements of stiffness matrix. (3) Free state: the node motion is not related to each other. The implement method of free state is that zero values are assigned to off-diagonal elements of stiffness matrix. Three main states are implemented through off-diagonal elements of stiffness matrix as shown in Table 1.

Table 1.

Properties of three main states of FEAO algorithm.

State	Off-diagonal elements	Nodes motion
Solid state	-	Limit
Diffusion state	+	Promote
Free state	0	No effect

Handling method of solving one-dimensional problem extends to multi-dimensional problem. By the principle of FEA element stiffness matrix, the diagonal element $k_{ii}$ ( $k_{jj}$ ) denotes the value of the force which is required to be exerted on node i(j), while node i(j) generates per unit displacement, with other nodes remaining 0 displacement. According to an optimization problem with n dimensions, stiffness equation (3) of element is described as follows: i, j——the node numbers; The off-diagonal element $k_{z, z + n}$ denotes value of the force which is required to be exerted on node i(j) while node j(i) generates per unit displacement in zth dimension, with other node remaining 0 displacement. Positive sign of $k_{z, z + n}$ , negative sign of $k_{z, z + n}$ and zero of $k_{z, z + n}$ correspond to diffusion state, solid state, and free state respectively.

K_{i j} = [\begin{matrix} \begin{matrix} 1 & 2 \dots & n \\ ↓ & ↓ & ↓ \end{matrix} & \begin{matrix} n + 1 & n + 2 \dots & 2 n \\ ↓ & ↓ & ↓ \end{matrix} \\ \begin{matrix} k_{i i} & \dots & 0 \\ ⋮ & ⋱ & ⋮ \\ 0 & \dots & k_{i i} \end{matrix} & \begin{matrix} k_{1, 1 + n} & \dots & 0 \\ ⋮ & ⋱ & ⋮ \\ 0 & \dots & k_{n, 2 n} \end{matrix} \\ \begin{matrix} k_{1 + n, 1} & \dots & 0 \\ ⋮ & ⋱ & ⋮ \\ 0 & \dots & k_{2 n, n} \end{matrix} & \begin{matrix} k_{j j} & \dots & 0 \\ ⋮ & ⋱ & ⋮ \\ 0 & \dots & k_{j j} \end{matrix} \end{matrix}]

(3)

An optimization based on FEA

The FEAO algorithm addresses the problem of locating the global minima of a (nonlinear) objective function f(x) defined on a bounded hyperspace; that is, determine:

$\min {f (X) : X \in$ Ω $\subset R^{n}}$ , $f : Ω \subset R^{n} \to R$ , $where Ω = {X \in Ω | x_{k}^{\min} \leq x_{k} \leq x_{k}^{\max}, k = 1, \dots, n}$ is the bounded region of feasible (allowable) solutions. n: dimensionality of the problem.

$x_{k}^{\max}$ : upper bound in each dimension. $x_{k}^{\min}$ : lower bound in each dimension. f(X): pointer to the function being minimized. $N_{pop}$ : population size.

FEAO framework

The FEAO algorithm comprises four main procedures: (a) Initialization;(b) Global stiffness matrix Calculation; (c) Nodes force calculation; (d) Displacement and position solution. FEAO’s algorithm framework (pseudocode) appears in Table 2.

Table 2.

FEAO framework.

BeginInitialize population: both displacement and position;Set parameters: population size

N_{pop}; \max_length

—— the maximum number of iterations for run termination;Iteration = 1;

While (termination criterion is not met)DoEvaluate all individuals using corresponding fitness function.Calculate the diagonal elements in stiffness element matrix of element established by every two individuals.Divide the node into three states, and calculate the off-diagonal elements in stiffness element matrix.Assembly element stiffness matrix to establish global stiffness matrix K.Calculate the node force vector F;solve stiffness equation

KS = F

to get node displacement S and update position X.Evaluate all individuals using corresponding fitness function;Iteration = t Iteration + 1;End Do (until termination criterion is met)End

Initializing population

In the Initialization step, $N_{pop}$ individuals are randomly generated in range $[x_{k}^{\min} x_{k}^{\max}]$ . The fitness (value of the objective function) for each individual is calculated using the function pointer f(X). $X_{best}$ —— the global best position at step t. $X_{i} (t) = (x_{i, 1} (t), x_{i, 2} (t) \dots x_{i, n} (t))$ ——position vector of node i at step t, $x_{i, k} (t)$ —— position component in k dimension. $U_{i} (t) = (u_{i, 1} (t), u_{i, 2} (t) \dots u_{i, n} (t))$ —— displacement vector of node i at step t, $u_{i, k} (t)$ ——displacement component of node i in kth dimension.

Calculating stiffness matrix

Calculating the diagonal elements in stiffness matrix

While any two node i and node j constitute a bar element, the stiffness coefficient $k_{ii}$ and $k_{jj}$ as the diagonal elements are calculated from the objective function f(x). Clearly it is possible to define the stiffness function using different methods. Adopting mass function in the APO algorithm,¹⁵ the stiffness function $k_{ii} = g (f (x_{i}))$ is defined as equation

\begin{matrix} G (f (x_{i})) = e^{\frac{f (x_{best}) - f (x_{i})}{f (x_{worst}) - f (x_{best}))}} \end{matrix}

(4)

According to the equation $k_{ii} = e^{\frac{f (x_{best}) - f (x_{i})}{f (x_{worst}) - f (x_{best}))}}$ , (i = 1.2…n) and $k_{jj} = e^{\frac{f (x_{best}) - f (x_{j})}{f (x_{worst}) - f (x_{best}))}}$ (j = n + 1…2n), $k_{ii}$ and $k_{jj}$ are assigned to the diagonal elements in stiffness matrix $K_{ij}$ as follow by equation (5).

K_{i j} = [\begin{matrix} \begin{matrix} 1 & \dots & n \\ ↓ & ↓ \end{matrix} & \begin{matrix} n + 1 & \dots & 2 n \\ ↓ & ↓ \end{matrix} \\ \begin{matrix} k_{i i} & \dots & 0 \\ ⋮ & ⋱ & ⋮ \\ 0 & \dots & k_{i i} \end{matrix} & \begin{matrix} 0 & \dots & 0 \\ ⋮ & ⋱ & ⋮ \\ 0 & \dots & 0 \end{matrix} \\ \begin{matrix} 0 & \dots & 0 \\ ⋮ & ⋱ & ⋮ \\ 0 & \dots & 0 \end{matrix} & \begin{matrix} k_{j j} & \dots & 0 \\ ⋮ & ⋱ & ⋮ \\ 0 & \dots & k_{j j} \end{matrix} \end{matrix}]

(5)

Calculating off-diagonal elements in stiffness matrix

According to the equation $cl a_{i} = \frac{f (x_{\max}) - f (x_{i})}{f (x_{\max}) - f (x_{\min})}$ , the method for node state division maps the fitness values into interval (0, 1) to be more reasonable in the division of states. When the condition $cl a_{i} >$ $standar d_{s}$ is satisfied, nodes may enter solid state in the dimension. When the condition $cl a_{i} \leq$ $standar d_{g}$ is satisfied, nodes may enter diffusion state in all dimensions. $Standar d_{s}$ is dividing standard for solid state. And $standar d_{g}$ is dividing standard for diffusion state. $Standar d_{s}$ and $standar d_{g}$ could be set according the function to be optimized.

A bar element which consist of node i and node j enters solid state when the conditions is simultaneously satisfied as follow: (1) $| x_{j, k} - x_{i, k} | \leq L$ , L—activated distance of solid state; (2) $cl a_{i} >$ $standar d_{s}$ or $cl a_{j} >$ $standar d_{s}$ , $standar d_{s} - - dividing standard of solid state$ ; (3) According to descending order of $cl a_{d} (d = 1, 2 \dots N_{pop})$ , $cl a_{i}$ or $cl a_{j}$ is within the top N( $1 < N \leq N_{pop}$ ). When the bar elements satisfy all above conditions, off-diagonal elements $k_{z, z + n} i = 1, 2 \dots n$ are assigned utilizing the equation (6). Thus, Off-diagonal element $k_{z, z + n} i = 1, 2 \dots n$ is set in element stiffness matrix $K_{ij}$ by equation (7).

\begin{matrix} k_{z, z + n} = - 1 \times \frac{T_{1} \times t}{\max_length}, z = 1, 2 \dots n \end{matrix}

(6)

$T_{1}$ ——the coefficient of solid force; $\max_length$ —— the maximum number of iterations for run termination; t——iterations.

\begin{matrix} K_{ij} = [\begin{matrix} \begin{matrix} k_{ii} & \dots & 0 \\ ⋮ & ⋱ & ⋮ \\ 0 & \dots & k_{ii} \end{matrix} & \begin{matrix} k_{1, 1 + n} & \dots & 0 \\ ⋮ & ⋱ & ⋮ \\ 0 & \dots & k_{n, 2 n} \end{matrix} \\ \begin{matrix} k_{1 + n, 1} & \dots & 0 \\ ⋮ & ⋱ & ⋮ \\ 0 & \dots & k_{2 n, n} \end{matrix} & \begin{matrix} k_{jj} & \dots & 0 \\ ⋮ & ⋱ & ⋮ \\ 0 & \dots & k_{jj} \end{matrix} \end{matrix}] \end{matrix}

(7)

A bar element which consist of node i and node j enters diffusion state when the conditions is simultaneously satisfied as follow: (1) $cl a_{i} \leq$ $standar d_{g}$ and $cl a_{j} \leq$ $standar d_{g}$ , $standar d_{g}$ — dividing standard of diffusion state (2) According to ascending order of $cl a_{d} (d = 1, 2 \dots N_{pop})$ , $cl a_{i}$ and $cl a_{j}$ are within the top N( $1 < N \leq N_{pop}$ ). When the bar elements satisfy all above conditions, off-diagonal elements $k_{z, z + n} (z = 1, 2 \dots n)$ are assigned utilizing the equation (8) as follow. Thus, Off-diagonal element $k_{z, z + n} i = 1, 2 \dots n$ is set in element stiffness matrix $K_{ij}$ utilizing equation (7).

\begin{matrix} k_{z, z + n} (z = 1, 2 \dots n) = (+ 1) \times \frac{T_{2} \times t}{\max_length} \end{matrix}

(8)

$T_{2}$ ——the coefficient of diffusion force; $\max_length$ —— the maximum number of iterations for run termination; t——iterations.

Assembling element stiffness matrix to establish global stiffness matrix K

K = [\begin{matrix} (N_{pop} - 1) k_{ii} & \dots & 0 & k_{1, m} & \dots & k_{1, n} \\ ⋮ & ⋱ & ⋮ & ⋮ & ⋱ & ⋮ \\ 0 & \dots & (N_{pop} - 1) k_{pp} & 0 & \dots & 0 \\ k_{q, 1} & \dots & 0 & (N_{pop} - 1) k_{qq} & \dots k_{q, r} & 0 \\ ⋮ & ⋱ & ⋮ & ⋮ & ⋱ & ⋮ \\ k_{n, 1} & \dots & k_{2 n, n} & 0 & k_{q, r} & (N_{pop} - 1) k_{jj} \end{matrix}] \begin{matrix} \leftarrow p \\ \leftarrow q \end{matrix}

(9)

Where p—free state, q—solid or diffusion state.

Because in pth row only the diagonal element $K_{pp} = (N_{pop} - 1) k_{pp}$ is a non-zero value, node p is in free state as equation (9) shows. Except the diagonal element $K_{pp} = (N_{pop} - 1)$ $k_{qq}$ in qth row, there are non-zero value in the off-diagonal elements, therefore node q is in solid or diffusion state.

Calculating the node force vector

Based on the results of the stiffness coefficient $k_{ii}$ and $k_{jj}$ , the node forces are calculated by an attraction- repulsion rule as follow.

\begin{matrix} P_{ij, k} = {\begin{matrix} G k_{ii} k_{jj} r_{ij . k} iff (X_{j}) < f (X_{i}) \\ G k_{ii} k_{jj} r_{ji . k} iff (X_{j}) \geq f (X_{i}) \end{matrix} \forall i \neq j \end{matrix}

(10)

$P_{ij, k}$ is the kth component of force which is exerted on node i by node j; $r_{ij . k} = x_{j, k} - x_{i, k}$ , $x_{i, k} and x_{j, k}$ are the kth-dimension coordinate. Equation (10) shows that node j attracts node i, and $P_{ij, k}$ is an attractive force if $f (X_{j}) < f (X_{i})$ ; or node j repels node i, $P_{ij, k}$ is a repulsive force if $f (X_{j}) \geq f (X_{i})$ .¹⁶

The kth component of the total force $P_{i, k}$ which is exerted on node i by all other nodes is obtained by summing over all other node forces as follow.

\begin{matrix} P_{i, k} = \sum_{\binom{j = 1,}{j \neq i}}^{N_{pop}} P_{ij, k} \end{matrix}

(11)

In this equation (11) of $P_{i, k}$ , random variable $α$ normally distributed on (0,1) is introduced.

\begin{matrix} P_{i, k} = α \sum_{\binom{j = 1,}{j \neq i}}^{N_{pop}} P_{ij, k} \end{matrix}

(12)

According to theory of FEA, the node force vector is calculated as follow.

P^{e} = {[P_{1, 1} P_{1, 2} \dots P_{1, n} \dots P_{2, 1} \dots P_{N_{pop}, 1} P_{N_{pop}, 2} \dots P_{N_{pop}, n}]}^{T}

Solving node displacement

Only the diagonal element of nodes in free state is a non-zero value. Thus, displacement of the free nodes could be obtained by the direct calculation, while displacement in solid or diffusion state are calculated through solving global stiffness equation.

Global stiffness equation (13) is as follow:

\begin{matrix} KS (t) = P^{e} \end{matrix}

(13)

\begin{matrix} S (t) = (s_{1, 1} (t), s_{1, 2} (t) \dots s_{1, n} (t) \dots s_{N_{pop}, 1} (t), s_{N_{pop}, 2} (t) \dots s_{N_{pop}, n} (t)) — — the displacement \\ vector calculated through stiffness equation \end{matrix}

$f (X_{best})$ is the objective function value at the position of the "best" node. According to the finite element method, the constraint condition is that the displacement is zero at the position of the “best” node. The constraint condition is imposed by crossing out the rows and columns which involve $X_{best}$ of global stiffness equation. $0 \leq ω < 1$ is a user-specified inertia weight that determines how easily the previous displacement can be changed, with larger values of $ω$ corresponding to greater displacement changes.

The final FEAO step is to update position, which refers to the displacement of nodes through the decision space using the previously computed total force to calculate a “displacement,” which then is used to update the nodes’ position. The displacement and coordinates of node i at the next iteration, are updated using equations (14) and (15), respectively:

\begin{matrix} u_{i, k} (t + 1) = ω u_{i, k} (t) + s_{i, k} (t) \end{matrix}

(14)

\begin{matrix} x_{i, k} (t + 1) = x_{i, k} (t) + u_{i, k} (t + 1) \end{matrix}

(15)

Parallel analysis

The FEAO algorithm adopts the parallel finite element method. Any element stiffness matrix is only dependent on themselves, and does not involve other element. All element stiffness matrix is simultaneously calculated, while the elements of matrix are implemented with parallelism.¹⁷ The global stiffness matrix is assembled by the element stiffness matrix in accordance with the principle of “matching the number of nodes”. In the phase of assembling matrix, the processes in which all elements that have no common nodes with each other are superimposed are handled concurrently. The global stiffness equation is solved through both direct and iterative method.¹⁸ There are mainly parallel direct methods: parallel Gauss-Jordan method, parallel Gauss method, parallel multifrontal method. Iterative method mainly refers to parallel preprocessing conjugate gradient method.^19–20

Convergence analysis

The paper analyzes the convergence of FEAO algorithm theoretically, adopting the method by which Zeng Jianchao et al proved the algorithm convergence of PSO²¹ and APO²² based on discrete-time linear system theory. Nodes which are in solid state and diffusion state not always exist in iterative process, but there are always nodes in free state. Convergence of the algorithm is proved by convergence analysis of nodes in free state. Only the diagonal element of nodes in free state is a non-zero value. Therefore, the global stiffness equations are simplified to the root equations of a degree. The solutions of the stiffness equation (13) simplified are as follow.

\begin{matrix} s_{i, k} (t) = α \sum_{\binom{j = 1}{j \neq i}}^{N_{pop}} \frac{P_{i, j}}{(N_{pop} - 1) k_{ii}} \end{matrix}

(16)

When an individual i chosen arbitrarily can prove the convergence, the result could extent to all the others. Each dimension of equations (10)–(15) is irrelevant to the other dimension. $X_{best}$ is assumed to be constant over some number of iterations. For reducing the convergence analysis, the paper drops the dimension index k in equations (14) and (15), the update equations become:

\begin{matrix} u_{i} (t + 1) = ω u_{i} (t) + α \sum_{\binom{j = 1}{j \neq i}}^{N_{pop}} \frac{P_{i, j}}{(N_{pop} - 1) k_{ii}} \end{matrix}

(17)

\begin{matrix} x_{i} (t + 1) = x_{i} (t) + u_{i} (t + 1) \end{matrix}

(18)

Define $N_{i} = {j | f (X_{j}) < f (X_{i}) \forall j \in S} and M_{i} = {j | f (X_{j}) \geq f (X_{i}) \forall j \in S}$ . Where, S is the set of all individuals. The displacement update equation (18) becomes:

\begin{array}{l} u_{i} (t + 1) = ω u_{i} (t) + α \sum_{\begin{matrix} j = 1 \\ j \neq i \end{matrix}}^{N_{p o p}} \frac{P_{i, k}}{(N_{p o p} - 1) k_{i i}} \\ = ω u_{i} (t) - α \sum_{j ϵ N_{i}} G k_{j j} (X_{i} (t) - X_{j} (t)) / (N_{p o p} - 1) \\ + α \sum_{j ϵ M_{i}} G k_{j j} (X_{i} (t) - X_{j} (t)) / (N_{p o p} - 1) \\ = ω u_{i} (t) + (α \sum_{j ϵ M_{i}} G k_{j j} - α \sum_{j ϵ N_{i}} G k_{j j}) X_{i} (t) / (N_{p o p} - 1) \\ - (α \sum_{j ϵ M_{i}} G k_{j j} - α \sum_{j ϵ N_{i}} G k_{j j}) X_{j} (t) / (N_{p o p} - 1) \end{array}

(19)

With the following definitions:

G_{N_{i}} = \frac{α \sum_{j ϵ N_{i}} G k_{jj}}{N_{pop} - 1}

(20)

\begin{matrix} G_{M_{i}} = \frac{α \sum_{j ϵ M_{i}} G k_{jj}}{N_{pop} - 1} \end{matrix}

(21)

\begin{matrix} G_{N M_{i}} = G_{N_{i}} - G_{M_{i}} \end{matrix}

(22)

\begin{matrix} F_{G_{i}} = \frac{α (\sum_{j ϵ N_{i}} G k_{jj} - \sum_{j ϵ M_{i}} G k_{jj}) X_{j} (t)}{N_{pop} - 1} \end{matrix}

(23)

Equation (19) becomes:

\begin{matrix} u_{i} (t + 1) = ω u_{i} (t) - G_{N M_{i}} X_{i} (t) + F_{G_{i}} \end{matrix}

(24)

Substituting equation (24) into equation (18) yields the following 2nd order difference equation.

\begin{matrix} X_{i} (t + 1) - (1 + ω - G_{N M_{i}}) X_{i} (t) + ω X (t - 1) = F_{G_{i}} \end{matrix}

(25)

Through property analysis of equation (25), the convergence of the RV sequence {E $[X_{i} (t)]$ } may be proved, E ${X_{i} (t)}$ ——the expectation of random variable $X_{i} (t)$ . Applying the expected value operator to equation (25), the equation (25) become:

\begin{matrix} E X_{i} (t + 1) - (1 + ω - E (G_{N M_{i}})) E X_{i} (t) + ω EX (t - 1) = E (F_{G_{i}}) \end{matrix}

(26)

where

\begin{matrix} φ_{i} = E (G_{N M_{i}}) = \frac{1}{2} G (\frac{\sum_{j ϵ N_{i}} k_{jj} - \sum_{j ϵ M_{i}} k_{jj}}{N_{pop} - 1}) \end{matrix}

(27)

Θ_{i} = E (F_{G_{i}}) = \frac{1}{2} G (\sum_{j ϵ N_{i}} k_{jj} E X_{j} (t) - \sum_{j ϵ M_{i}} k_{jj} E X_{j} (t)) \frac{1}{N_{pop} - 1}

(28)

The characteristic equation of equation (25) is:

\begin{matrix} λ^{2} - (1 + ω - φ_{i}) λ + ω = 0 \end{matrix}

(29)

Theorem 1. If and only if 0 ≤ $ω$ <1 and 0 $< φ_{i} < 2 (1 + ω)$ converges to $X_{best}$ .

Proof. The condition of { $E X_{i} (t)$ } convergence is that the absolute values of both eigenvalues are less than 1.

\begin{matrix} | \frac{1 + ω - φ_{i} \pm \sqrt{{(1 + ω - φ_{i})}^{2} - 4 ω}}{2} | < 1 \end{matrix}

(30)

There are three cases that depend on the value of the discriminant ${(1 + ω - φ_{i})}^{2} - 4 ω$ as follows.

(1) ${(1 + ω - φ_{i})}^{2} = 4 ω$ , The eigenvalues are real and equal, (2) ${(1 + ω - φ_{i})}^{2} < 4 ω$ , both eigenvalues are real, but different. (3) ${(1 + ω - φ_{i})}^{2} > 4 ω$ , the two eigenvalues are complex conjugates with equal magnitudes.

By summarizing three cases, the convergence condition for the sequence of random variables {E $X_{i} (t)$ } is obtained as follow.

0 \leq ω < 1 and 0 < φ_{i} < 2 (1 + ω)

(31)

When {E $X_{i} (t)$ } is convergent, with the assumption that limit of convergence is $X^{a}$ , substituting $lim_{t \to \infty} E X_{i} (t) = X^{a}$ into equation (26) yields

\begin{matrix} φ_{i} X^{a} = θ_{i} \end{matrix}

(32)

Substituting equations (27) and (28) into equation (32) yields

\frac{1}{2} G \frac{(\sum_{j ϵ N_{i}} k_{jj} - \sum_{j ϵ M_{i}} k_{jj}) X^{a}}{N_{pop} - 1} = \frac{1}{2} G \frac{\sum_{j ϵ N_{i}} k_{jj} E X_{j} (t) - \sum_{j ϵ M_{i}} k_{jj} X_{j} (t)}{N_{pop} - 1}

which may be written

\begin{matrix} \Rightarrow \sum_{j ϵ N_{i}} k_{jj} X^{a} - \sum_{j ϵ M_{i}} k_{jj} X^{a} - \sum_{j ϵ N_{i}} k_{jj} E X_{j} (t) + \sum_{j ϵ M_{i}} k_{jj} E X_{j} (t)) = 0 \\ \Rightarrow \sum_{j ϵ N_{i}} k_{jj} (X^{a} - E X_{j} (t)) - \sum_{j ϵ M_{i}} k_{jj} (X^{a} - E X_{j} (t))) = 0 \\ \Rightarrow k_{best} (X^{a} - E X_{j} (t)) + \sum_{\begin{matrix} j ϵ N_{i} \\ j \neq best \end{matrix}} k_{jj} (X^{a} - E X_{j} (t)) - \sum_{j ϵ M_{i}} k_{jj} (X^{a} - E X_{j} (t)) = 0 \end{matrix}

(33)

By inspection, equation (33) is satisfied if and only $E X_{j} (t) = X^{a} = X_{best}$ . Substituting equation (27) into equation (31) yields

0 \leq ω < 1 and 0 < G (\sum_{j ϵ N_{i}} k_{jj} - \sum_{j ϵ M_{i}} k_{jj}) < 4 (1 + ω) (N_{pop} - 1)

(34)

According to equation (34), the algorithm converges by setting an appropriately small enough value for G.

Numerical experiments

In this section, FEAO’s performance is compared with PSO’s through the six recognized benchmark functions listed in Table 3. The Table 3 includes the search domain, the known minimum function value.

Table 3.

Parameters for Test Function.

Function	n	$N_{pop}$	[ $x^{\min}, x^{\max}$ ]
Tablet	30	20	[–100,100]
Quadric	30	20	[–100,100]
Rosenbrock	30	20	[–50,50]
Griewank	30	20	[–300,300]
Rastrigin	30	20	[–5.12,5.12]
Schaffer’s f7	30	20	[–100,100]

FEAO algorithm runs were made with the following empirically determined parameters: inertia weight $ω = 0.9 - \frac{t}{MAXITER} \times 0.5$ , where t is iteration number and MAXITER is the maximum number of iterations for run termination, MAXITER = 6000; constant G = 0.29, problem dimensionality n = 30; population size $N_{pop}$ = 20; activated distance of solid state L = 1.0e-05;Solid coefficient $T_{1}$ = 5, diffusion coefficient $T_{2}$ = 3; diffusion threshold $standar d_{g} = 1$ ;solid threshold $standar d_{s} = 0.7$ ; $C_{1}$ ——the terminal iteration of solid state; $C_{2}$ ——the start iteration of diffusion state; N——the maximum number of nodes in solid or diffusion state.

Each numerical experiment was repeated 50 times. The recorded performance data were the best function value (BEST), the average best function value (Mean), the max function value (MAX), and its standard deviation (STD). The PSO test results were selected from the references.²³ The performance of FEAO algorithm is compared with PSO for five types which include free state, diffusion state, solid state, entirely solid state, synthesis state. Free state——in whole iterative process (t< MAXITER) all nodes are in free state; Diffusion state ——later in the iteration process (t> $C_{2}$ ) part of nodes are in diffusion state; Solid state ——Early in the iteration process (t< $C_{1}$ ) part of nodes are in solid state; Entirely solid state——in iterative process (t< MAXITER) part of nodes are in solid state; Synthesis state——Early in the iteration process (t< $C_{1}$ ) part of nodes are in solid state and later in the iteration process (t> $C_{2}$ ) part of nodes are in diffusion state.

Table 4 shows the comparative results. The results of testing Tablet function illustrate that except entirely solid state the other diffusion state, solid state and synthesis state of the FEAO has higher performance than the free state, and synthesis state obtains the best performance outperforming PSO. For Rosenbrock function, the numeral experiment indicates that the other four states have the advantage over free state, and synthesis state is the best exceeding PSO. For Quadric function, through simulations, compared with PSO and free state, the other four states enhance greatly the performance. For Griewank function, the performance of free state, diffusion state, entirely solid state, and synthesis state are similar to PSO, and solid state is the best exceeding the others. By the Rastrigin and Schaffer’s f7 test functions, the searching ability of the solid state and entirely solid state is more efficient, and the performance in all states to test Schaffer’s f7 are better than PSO.

Table 4.

Shows the comparative results.

Function	State	Best	Paramaters	Mean	MAX	STD
Tablet	PSO	8.2e-22		1.5e-19	5.6e-18	7.8e-19
	free	4.23e-13	$C_{1}$ = 3000; $C_{2}$ = 4000; $T_{1}$ = 5; $T_{2}$ = 3;L = 0.00001; $standar d_{s}$ = 0.7;N = 10; $standar d_{g}$ = 1;	1.04e-11	1.35e-09	3.1e-10
	diffusion	3.88e-22		1.60e-20	6.04e-18	9.8647e-19
	solid	4.13e-19		5.71e-16	9.45e-15	3.3219e-15
	entirely solid	9.00e-06		0.00125	0.0700	0.019
	synthesis	6.66e-25		7.83e-21	5.88e-20	8.64e-21
Rosenbrock	PSO	4.873		93.71	1.0e + 03	1.5e + 02
	free	104.56	$C_{1} = 3000$ ; $C_{2}$ = 4000;The remaining parameters are the same as above;	658.4	5.56e + 3	835.5
	diffusion	56.9		316.8	6.21e + 03	9.21e + 2
	solid	65.5		116.5	1508.5	375.5
	entirely solid	15.0		359.0	750	250.2
	synthesis	14.5		25.74	120.5	24.89
Quadric	PSO	0.772		1.0e + 2	4.8e + 02	1.0e + 2
	free	6.35e-13	$C_{1} = 3000$ ; $C_{2} = 4000$ ;The remaining parameters are the same as above;	4.2e-06	1.82	0.32
	diffusion	3.36e-26		7.96e-23	2.78e-20	4.9e-21
	solid	8.56e-22		1.8e-19	8.79e-9	2.01e-9
	entirely solid	6.12e-15		1.12e-12	6.52e-9	9.8e-10
	synthesis	4.56e-33		2.66e-31	5.6e-22	1.2e-22
Griewank	PSO	6.7e-11		2.6e-02	0.108	2.8e-02
	free	2.56e-7	$C_{1} = 3000$ ; $C_{2} = 4000$ ;The remaining parameters are the same as above;	0.012	0.095	3.2e-02
	diffusion	9.84e-10		0.012	0. 13	4e-02
	solid	2.45e-18		2.22e-16	8.75e-13	4.34e-13
	entirely solid	9.0e-8		0.0073	0.02	0.005
	synthesis	5.6e-10		0.063	0.25	0.07
Rastrigin	PSO	30.84		55.18	81.59	12.31
	free	21.8	The same as above;	41.79	107	15.9
	diffusion	23.7		37.81	98.5	12.3
	solid	15.5		22.88	89.6	11.9
	entirely solid	10.4		14.10	55.8	9.8
	synthesis	29.8		42.78	78.6	6.32
Schaffer’s f7	PSO	2.4e + 02		2.9e + 02	3.5e + 02	30.91
	free	0.26	The same as above;	0.42	12.4	3.2
	diffusion	0.37		1.29	9.08	2.6
	solid	0.20		0.27	4.95	1.32
	entirely solid	0.20		0.28	3.77	0.97
	synthesis	0.78		1.59	17.9	4.68

The paper compared the dynamic performance of five states which are sampled at an interval of 20 iterations as shown in Figures 3 –8. The results manifest that the diffusion state, solid state and synthesis state have advantages over free state. Because the force is proportional to the distance of nodes, at the end of the iteration the distance of two nodes is too close to generate large enough force for escaping from local minima. Therefore, the dynamic curve of the free state levels off as show. Each state has the capabilities and limitations. Solid state is well-suited for dealing with the Griewank, Rastrigin, Schaffer’s f7 functions, while synthesis and diffusion states are fitter for solving the problem of Tablet, Rosenbrock, Quadric. Solid state can improve the diversity of the population, and diffusion state enhances node force to improve search performance.

Figure 3.

Dynamic performance for Tablet.

Figure 4.

Dynamic performance for Rosenbrock.

Figure 5.

Dynamic performance for Quadric.

Figure 6.

Dynamic performance for Griewank.

Figure 7.

Dynamic performance for Rastrigin.

Figure 8.

Dynamic performance for Schaffer’s f7.

Conclusions and future work

Based on Finite Element Analysis, a novel optimization is proposed. The FEA approach that can simulate multi-states of matter is adopted to balance the global search ability and local search ability in the novel optimization algorithms. The condition for FEAO convergence is deduced through analyzing FEAO algorithm based on discrete-time linear system theory. The simulation results of six typical optimization test functions show that FEAO algorithm is effective. The paper analyzes parallelism at main steps by theoretical analysis. Thus, a further study of FEAO algorithm is needed to focus on realization of parallel/distributed computation.

Footnotes

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This project is supported by National Natural Science Foundation of China (Grant No. 51875381), National Natural Youth Science Foundation of China (Grant No. 51905369), and administration committee of Yangquan development zone (No.17-015).

ORCID iD

Youshan Gao

References

Holland

Adaptation in natural and artificial systems. Cambridge, MA: MIT press, 1992.

Eberhart

Kennedy

New optimizer using particle swarm theory. In: MHS’95 Proceedings of the sixth international symposium on micro machine and human science, Nagoya, Japan, 4–6 October 1995, pp. 39–43. Piscataway, NJ: IEEE.

Eberhart

Shi

Tracking and optimizing dynamic systems with particle swarms. In: Proceedings of the 2001 congress on evolutionary computation, Seoul, South Korea, 27–30 May 2001, pp. 94–100. Piscataway, NJ: IEEE.

Shen

Wang

, et al. A dynamic adaptive particle swarm optimization for knapsack problem. In: Proceedings of Wcica 2006: sixth world congress on intelligent control and automation, Dalian, China, 21–23 June 2006, pp. 3183–3187. IEEE.

Formato

. Central force optimization: a new nature inspired computational framework for multidimensional search and optimization. In: Krasnogor

Nicosia

Pavone

, et al. (eds) Nature inspired cooperative strategies for optimization. Berlin: Springer, 2008, 221–238.

Formato

Improved CFO algorithm for antenna optimization. Prog Electromagn Res 2010; 19: 405–425.

Formato

Central force optimization: a new deterministic gradient-like optimization metaheuristic. OPSEARCH 2009; 46(1): 25–51.

Birbil

Fang

RS.

An electromagnetism-like mechanism for global optimization. J Global Optim 2003; 25(3): 263–282.

Xie

Zeng

Cui

ZH.

On mass effects to artificial physics optimization algorithm for global optimization problems. Int J Innov Comput Applic 2010; 2(2): 69–76.

10.

Bathe

KJ.

Finite element procedures. Watertown, MA: Klaus-Jürgen Bathe, 2014.

11.

Noor

AK.

Parallel processing in finite element structural analysis. Eng Comput 1988; 3: 225–241.

12.

Kohei

Takuma

Tsuyoshi

, et al. Acceleration of element by element kernel in unstructured implicit low-order finite-element earthquake simulation using open ACC on Pascal GPUs. In: Third workshop on accelerator programming using directives, Salt Lake City, UT, USA, 14 November 2016, pp. 1–12. IEEE.

13.

Chung

Efendiev

Lee

CS.

Mixed generalized multiscale finite element methods and applications. Multiscale Model Sim 2015; 13(1): 338–366.

14.

Zeng

Fundamentals of finite element analysis. Beijing: Tsinghua University Press, 2008.

15.

Xie

Zeng

JC.

The performance analysis of artificial physics optimization algorithm driven by different virtual forces. ICIC-EL 2010; 4(1): 239–244.

16.

Maenhout

Vanhoucke

An electromagnetic meta-heuristic for the nurse scheduling problem. J Heuristics 2007; 13(4): 359–385.

17.

Yan

Tang

, et al. Parallel realization of element by element analysis of eddy current field based on graphic processing unit. Appl Comput Electrom Soc J 2018; 33(2): 168–171.

18.

Hsieh

KT.

Parallelization of MAP3D based on element-by-element Jacobi preconditioned conjugate gradient method. IEEE Trans Magn 2003; 39(1): 139–147.

19.

Meijerink

Vander Vorst

HA.

An iterative solution method for linear systems of which the coefficient matrix is a symmetric M-matrix. Math Comput 1977; 137(31): 148–162.

20.

Parallel computing of electromagnetic field based on domain decomposition method. In: Preprints of the 5th international conference on electric utility deregulation and restructuring and power technologies. Changsha, China, 26–29 November 2015, pp. 1583–1586. IEEE.

21.

Cui

Zeng

JC.

A new model of particle swarm optimization and its theoretical analysis. J Comput Res Develop 2006; 43(1): 96–100.

22.

Xie

Zeng

Formato

RA.

Convergence analysis and performance of the extended artificial physics optimization algorithm. Appl Math Comput 2011; 218(8): 4000–4011.

23.

Wang

Qing

, et al. An improved particle swarm optimization based on self-adaptive escapevelocity. J Softw 2005; 16(12): 2036-2044.

Design and performance analysis of a new optimization algorithm based on Finite Element Analysis

Abstract

Keywords

Introduction

Finite element analysis

A relationship of mapping between FEA and a population-based optimization algorithm

An optimization based on FEA

FEAO framework

Initializing population

Calculating stiffness matrix

Calculating the diagonal elements in stiffness matrix

Calculating off-diagonal elements in stiffness matrix

Assembling element stiffness matrix to establish global stiffness matrix K

Calculating the node force vector

Solving node displacement

Parallel analysis

Convergence analysis

Numerical experiments

Conclusions and future work

Footnotes

Declaration of conflicting interests

Funding

ORCID iD

References