Self-adaptive Potential-based Stopping Criteria for Particle Swarm Optimization

TL;DR

This paper introduces a self-adaptive stopping criterion for Particle Swarm Optimization based on counting forced moves, enabling the algorithm to terminate efficiently near local optima.

Contribution

It proposes a novel potential-based stopping criterion for PSO that self-adapts by monitoring forced moves, improving convergence speed.

Findings

Stopping criteria effectively detect convergence near optima.

The method finds good solutions faster than traditional criteria.

Theoretical analysis confirms behavior near local optima.

Abstract

We study the variant of Particle Swarm Optimization (PSO) that applies random velocities in a dimension instead of the regular velocity update equations as soon as the so-called potential of the swarm falls below a certain bound in this dimension, arbitrarily set by the user. In this case, the swarm performs a forced move. In this paper, we are interested in how, by counting the forced moves, the swarm can decide for itself to stop its movement because it is improbable to find better solution candidates as it already has found. We formally prove that when the swarm is close to a (local) optimum, it behaves like a blind-searching cloud, and that the frequency of forced moves exceeds a certain, objective function-independent value. Based on this observation, we define stopping criteria and evaluate them experimentally showing that good solution candidates can be found much faster than applying other criteria.

Tables

Table 1. Table 1 : Experimental results for ( γ , μ ) 𝛾 𝜇 (\gamma,\mu) - full stop . iter term subscript iter term \text{iter}_{\rm term} denotes the number of iterations (measured every μ 𝜇 \mu th iteration) until ( γ , μ ) 𝛾 𝜇 (\gamma,\mu) - full stop terminates the execution of f-PSO. G t subscript 𝐺 𝑡 G_{t} denotes the global attractor after t 𝑡 t iterations. The final two columns present the median and (with ± plus-or-minus \pm ) the standard deviation of the gradient at the global attractor over 500 500 500 runs of the f-PSO with and without ( γ , μ ) 𝛾 𝜇 (\gamma,\mu) - full stop . In the tests, N = 5 𝑁 5 N=5 , D = 15 𝐷 15 D=15 , μ = 50 000 𝜇 50000 \mu=50\,000 , σ stag ​ ( N , D , μ ) = 318 350 subscript 𝜎 stag 𝑁 𝐷 𝜇 318350 \sigma_{\rm stag}(N,D,\mu)=318\,350 , γ = 1350 𝛾 1350 \gamma=1350 , and δ = 10 − 7 𝛿 superscript 10 7 \delta=10^{-7} . The less the gradient, the better G t subscript 𝐺 𝑡 G_{t} .

Function $f$	${\text{iter}}_{\mathrm{term}}$	$\Vert \nabla f\left(G_{{\mathrm{iter}}_{\mathrm{term}}}\right)\Vert$	$\Vert \nabla f\left(G_{\mathrm{15\hspace{0.17em}000\hspace{0.17em}000}}\right)\Vert$
Sphere	$\mathrm{100\hspace{0.17em}000}\pm 0$	$6.65\times {10}^{-8}\pm 6.93\times {10}^{-9}$	$3.60\times {10}^{-8}\pm 3.83\times {10}^{-9}$
H. C. Elliptic	$\mathrm{300\hspace{0.17em}000}\pm \mathrm{154\hspace{0.17em}324}$	$2.22\times {10}^{-5}\pm 1.93\times {10}^{-5}$	$1.17\times {10}^{-5}\pm 1.09\times {10}^{-5}$
Schwefel	$\mathrm{150\hspace{0.17em}000}\pm \mathrm{110\hspace{0.17em}707}$	$1.94\times {10}^{-7}\pm 1.20\times {10}^{-7}$	$1.07\times {10}^{-7}\pm 5.96\times {10}^{-8}$
Rastrigin	$\mathrm{100\hspace{0.17em}000}\pm 5000$	$1.34\times {10}^{-5}\pm 1.61\times {10}^{-6}$	$7.82\times {10}^{-6}\pm 1.15\times {10}^{-6}$
Rosenbrock	$\mathrm{850\hspace{0.17em}000}\pm \mathrm{490\hspace{0.17em}538}$	$3.09\times {10}^{-5}\pm 6.45\times {10}^{-6}$	$2.04\times {10}^{-5}\pm 3.80\times {10}^{-6}$

Table 2. Table 2 : Experimental results for ( κ , γ , μ ) 𝜅 𝛾 𝜇 (\kappa,\gamma,\mu) - partial stop (for notions, see Table  1 ). Median and (with ± plus-or-minus \pm ) the standard deviation over 500 runs of the stopped f-PSO with N = 5 𝑁 5 N=5 , D = 15 𝐷 15 D=15 , μ = 50 000 𝜇 50000 \mu=50\,000 , σ stag ​ ( N , D , μ ) = 318 350 subscript 𝜎 stag 𝑁 𝐷 𝜇 318350 \sigma_{\rm stag}(N,D,\mu)=318\,350 , γ = 1350 𝛾 1350 \gamma=1350 , and δ = 10 − 7 𝛿 superscript 10 7 \delta=10^{-7} , and κ = 2 𝜅 2 \kappa=2 and κ = 8 𝜅 8 \kappa=8 , resp., against an ‘unstopped’ f-PSO terminated after 15 000 000 15000000 15\,000\,000 iterations.

Function $f$	${\text{iter}}_{\kappa =2}$	$\Vert \nabla f\left(G_{{\text{iter}}_{\kappa =2}}\right)\Vert$	${\text{iter}}_{\kappa =8}$	$\Vert \nabla f\left(G_{{\text{iter}}_{\kappa =8}}\right)\Vert$	$\Vert \nabla f\left(G_{\mathrm{15\hspace{0.17em}000\hspace{0.17em}000}}\right)\Vert$
Sphere	$\mathrm{50\hspace{0.17em}000}$ $\pm 0$	$7.24\times {10}^{-8}$ $\pm 7.50\times {10}^{-9}$	$\mathrm{50\hspace{0.17em}000}$ $\pm 0$	$7.27\times {10}^{-8}$ $\pm 8.17\times {10}^{-9}$	$3.60\times {10}^{-8}$ $\pm 3.83\times {10}^{-9}$
H. C. Elliptic	$\mathrm{50\hspace{0.17em}000}$ $\pm 0$	$2.52\times {10}^{-5}$ $\pm 2.22\times {10}^{-5}$	$\mathrm{50\hspace{0.17em}000}$ $\pm 0$	$2.44\times {10}^{-5}$ $\pm 2.24\times {10}^{-5}$	$1.17\times {10}^{-5}$ $\pm 1.09\times {10}^{-5}$
Schwefel	$\mathrm{50\hspace{0.17em}000}$ $\pm 0$	$2.22\times {10}^{-7}$ $\pm 1.24\times {10}^{-7}$	$\mathrm{50\hspace{0.17em}000}$ $\pm 0$	$2.14\times {10}^{-7}$ $\pm 1.31\times {10}^{-7}$	$1.07\times {10}^{-7}$ $\pm 5.96\times {10}^{-8}$
Rastrigin	$\mathrm{50\hspace{0.17em}000}$ $\pm 0$	$1.45\times {10}^{-5}$ $\pm 1.65\times {10}^{-6}$	$\mathrm{50\hspace{0.17em}000}$ $\pm 0$	$1.46\times {10}^{-5}$ $\pm 1.72\times {10}^{-6}$	$7.82\times {10}^{-6}$ $\pm 1.15\times {10}^{-6}$
Rosenbrock	$\mathrm{100\hspace{0.17em}000}$ $\pm \mathrm{730\hspace{0.17em}696}$	$8.54\times {10}^{-5}$ $\pm 3.78\times {10}^{-5}$	$\mathrm{150\hspace{0.17em}000}$ $\pm \mathrm{469\hspace{0.17em}614}$	$7.02\times {10}^{-5}$ $\pm 1.80\times {10}^{-5}$	$2.04\times {10}^{-5}$ $\pm 3.80\times {10}^{-6}$

Table 3. Table 3 : Experimental results for ( κ , γ , μ ) 𝜅 𝛾 𝜇 (\kappa,\gamma,\mu) - partial stop (for notions, see Table  1 ). Geometric Mean over 500 runs of the stopped f-PSO with N = 5 𝑁 5 N=5 , D = 15 𝐷 15 D=15 , μ = 50 000 𝜇 50000 \mu=50\,000 , σ stag ​ ( N , D , μ ) = 318 350 subscript 𝜎 stag 𝑁 𝐷 𝜇 318350 \sigma_{\rm stag}(N,D,\mu)=318\,350 , γ = 1350 𝛾 1350 \gamma=1350 , and δ = 10 − 7 𝛿 superscript 10 7 \delta=10^{-7} , and κ = 2 𝜅 2 \kappa=2 and κ = 8 𝜅 8 \kappa=8 , resp., against an ‘unstopped’ f-PSO terminated after 15 000 000 15000000 15\,000\,000 iterations.

Function $f$	${\text{iter}}_{\kappa =2}$	$\Vert \nabla f\left(G_{{\text{iter}}_{\kappa =2}}\right)\Vert$	${\text{iter}}_{\kappa =8}$	$\Vert \nabla f\left(G_{{\text{iter}}_{\kappa =8}}\right)\Vert$	$\Vert \nabla f\left(G_{\mathrm{15\hspace{0.17em}000\hspace{0.17em}000}}\right)\Vert$
Sphere	$\mathrm{50\hspace{0.17em}000}$	$7.15\times {10}^{-8}$	$\mathrm{50\hspace{0.17em}000}$	$7.12\times {10}^{-8}$	$3.60\times {10}^{-8}$
H. C. Elliptic	$\mathrm{50\hspace{0.17em}000}$	$1.95\times {10}^{-5}$	$\mathrm{50\hspace{0.17em}000}$	$1.95\times {10}^{-5}$	$1.17\times {10}^{-5}$
Schwefel	$\mathrm{50\hspace{0.17em}000}$	$1.04\times {10}^{-7}$	$\mathrm{50\hspace{0.17em}000}$	$1.07\times {10}^{-7}$	$1.07\times {10}^{-7}$
Rastrigin	$\mathrm{50\hspace{0.17em}000}$	$1.43\times {10}^{-5}$	$\mathrm{50\hspace{0.17em}000}$	$1.44\times {10}^{-5}$	$7.82\times {10}^{-6}$
Rosenbrock	$\mathrm{153\hspace{0.17em}672}$	$8.80\times {10}^{-5}$	$\mathrm{176\hspace{0.17em}989}$	$6.91\times {10}^{-5}$	$2.04\times {10}^{-5}$

Table 4. Table 4 : Experimental results for ( κ , γ , μ ) 𝜅 𝛾 𝜇 (\kappa,\gamma,\mu) - partial stop with reduced interval length (for notions, see Table  1 ). Median and (with ± plus-or-minus \pm ) the standard deviation over 500 runs of the stopped f-PSO with N = 5 𝑁 5 N=5 , D = 15 𝐷 15 D=15 , 𝝁 = 5 000 𝝁 5000 \boldsymbol{\mu=5\,000} , 𝝈 stag ​ ( 𝑵 , 𝑫 , 𝝁 ) = 31 835 subscript 𝝈 stag 𝑵 𝑫 𝝁 31835 \boldsymbol{\sigma}_{\textbf{stag}}\boldsymbol{(N,D,\mu)=31\,835} , γ = 1350 𝛾 1350 \gamma=1350 , and δ = 10 − 7 𝛿 superscript 10 7 \delta=10^{-7} , and κ = 2 𝜅 2 \kappa=2 and κ = 8 𝜅 8 \kappa=8 , resp., against an ‘unstopped’ f-PSO terminated after 15 000 000 15000000 15\,000\,000 iterations.

Function $f$	${\text{iter}}_{\kappa =2}$	$\Vert \nabla f\left(G_{{\text{iter}}_{\kappa =2}}\right)\Vert$	${\text{iter}}_{\kappa =8}$	$\Vert \nabla f\left(G_{{\text{iter}}_{\kappa =8}}\right)\Vert$	$\Vert \nabla f\left(G_{\mathrm{15\hspace{0.17em}000\hspace{0.17em}000}}\right)\Vert$
Sphere	$\mathrm{5\hspace{0.17em}000}$ $\pm 0$	$1.10\times {10}^{-7}$ $\pm 1.37\times {10}^{-8}$	$\mathrm{10\hspace{0.17em}000}$ $\pm 2200$	$9.52\times {10}^{-8}$ $\pm 1.22\times {10}^{-8}$	$3.63\times {10}^{-8}$ $\pm 3.73\times {10}^{-9}$
H. C. Elliptic	$\mathrm{5\hspace{0.17em}000}$ $\pm 1310$	$1.11\times {10}^{-4}$ $\pm 1.54\times {10}^{-4}$	$\mathrm{10\hspace{0.17em}000}$ $\pm 0$	$3.93\times {10}^{-5}$ $\pm 3.73\times {10}^{-5}$	$1.28\times {10}^{-5}$ $\pm 1.00\times {10}^{-5}$
Schwefel	$\mathrm{15\hspace{0.17em}000}$ $\pm 2502$	$3.87\times {10}^{-7}$ $\pm 2.18\times {10}^{-7}$	$\mathrm{15\hspace{0.17em}000}$ $\pm 1034$	$3.12\times {10}^{-7}$ $\pm 1.82\times {10}^{-7}$	$1.01\times {10}^{-7}$ $\pm 6.02\times {10}^{-8}$
Rastrigin	$\mathrm{5\hspace{0.17em}000}$ $\pm 0$	$2.12\times {10}^{-5}$ $\pm 2.69\times {10}^{-6}$	$\mathrm{5\hspace{0.17em}000}$ $\pm 2449$	$1.97\times {10}^{-5}$ $\pm 2.64\times {10}^{-6}$	$7.82\times {10}^{-6}$ $\pm 1.13\times {10}^{-6}$
Rosenbrock	$\mathrm{65\hspace{0.17em}000}$ $\pm \mathrm{442\hspace{0.17em}699}$	$1.40\times {10}^{-3}$ $\pm 1.04\times {10}^{-3}$	$\mathrm{80\hspace{0.17em}000}$ $\pm \mathrm{381\hspace{0.17em}917}$	$1.74\times {10}^{-4}$ $\pm 5.34\times {10}^{-5}$	$2.04\times {10}^{-5}$ $\pm 4.03\times {10}^{-6}$

Table 5. Table 5 : Experimental results for ( κ , γ , μ ) 𝜅 𝛾 𝜇 (\kappa,\gamma,\mu) - partial stop with reduced interval length (for notions, see Table  1 ). Geometric Mean over 500 runs of the stopped f-PSO with N = 5 𝑁 5 N=5 , D = 15 𝐷 15 D=15 , 𝝁 = 5 000 𝝁 5000 \boldsymbol{\mu=5\,000} , 𝝈 stag ​ ( 𝑵 , 𝑫 , 𝝁 ) = 31 835 subscript 𝝈 stag 𝑵 𝑫 𝝁 31835 \boldsymbol{\sigma}_{\textbf{stag}}\boldsymbol{(N,D,\mu)=31\,835} , γ = 1350 𝛾 1350 \gamma=1350 , and δ = 10 − 7 𝛿 superscript 10 7 \delta=10^{-7} , and κ = 2 𝜅 2 \kappa=2 and κ = 8 𝜅 8 \kappa=8 , resp., against an ‘unstopped’ f-PSO terminated after 15 000 000 15000000 15\,000\,000 iterations.

Function $f$	${\text{iter}}_{\kappa =2}$	$\Vert \nabla f\left(G_{{\text{iter}}_{\kappa =2}}\right)\Vert$	${\text{iter}}_{\kappa =8}$	$\Vert \nabla f\left(G_{{\text{iter}}_{\kappa =8}}\right)\Vert$	$\Vert \nabla f\left(G_{\mathrm{15\hspace{0.17em}000\hspace{0.17em}000}}\right)\Vert$
Sphere	$\mathrm{5\hspace{0.17em}000}$	$1.09\times {10}^{-7}$	$\mathrm{8\hspace{0.17em}337}$	$9.41\times {10}^{-8}$	$3.58\times {10}^{-8}$
H. C. Elliptic	$\mathrm{5\hspace{0.17em}262}$	$9.68\times {10}^{-5}$	$\mathrm{10\hspace{0.17em}000}$	$3.15\times {10}^{-5}$	$9.62\times {10}^{-6}$
Schwefel	$\mathrm{12\hspace{0.17em}276}$	$3.84\times {10}^{-7}$	$\mathrm{15\hspace{0.17em}150}$	$3.15\times {10}^{-7}$	$1.04\times {10}^{-7}$
Rastrigin	$\mathrm{5\hspace{0.17em}000}$	$2.09\times {10}^{-5}$	$\mathrm{6\hspace{0.17em}587}$	$1.96\times {10}^{-5}$	$7.75\times {10}^{-6}$
Rosenbrock	$\mathrm{92\hspace{0.17em}560}$	$1.41\times {10}^{-3}$	$\mathrm{111\hspace{0.17em}448}$	$1.80\times {10}^{-4}$	$2.01\times {10}^{-5}$

Equations

$$V^n$$

$$X^n$$

$${\Phi }_d:=\sum _{n=1}^N(\underset{=:{\varphi }_d^n}{\underbrace{|V_d^n|+|G_d-X_d^n|}})=\sum _{n=1}^N{\varphi }_d^n$$

$$V_d^n$$

$$X^n$$

$$\mathrm{P}\left[F_{i,2}=1|F_{i,1}=1\right]=\mathrm{P}\left[{\varphi }^{1,i-1}<\delta |F_{i,1}=1\right]$$

$$\begin{array}{rlrl}=& \mathrm{P}\left[{\varphi }^{1,i+1}<\delta |\Delta \ge 0\wedge F_{i,1}=1\right]& & \cdot \mathrm{P}\left[\Delta \ge 0|F_{i,1}=1\right]\\ +& \mathrm{P}\left[{\varphi }^{1,i+1}<\delta |\Delta <0\wedge F_{i,1}=1\right]& & \cdot \mathrm{P}\left[\Delta <0|F_{i,1}=1\right]\end{array}$$

$$\begin{array}{rlrl}=& \mathrm{P}\left[{\varphi }^{1,i+1}<\delta |(V^{1,i+1}\ge \Delta \ge 0)\wedge F_{i,1}=1\right]& & \cdot \mathrm{P}\left[V^{1,i+1}\ge \Delta \ge 0|F_{i,1}=1\right]\\ +& \mathrm{P}\left[{\varphi }^{1,i+1}<\delta |(\Delta \ge V^{1,i+1}\ge 0)\wedge F_{i,1}=1\right]& & \cdot \mathrm{P}\left[\Delta \ge V^{1,i+1}\ge 0|F_{i,1}=1\right]\\ +& \mathrm{P}\left[{\varphi }^{1,i+1}<\delta |(\Delta \ge 0\ge V^{1,i+1})\wedge F_{i,1}=1\right]& & \cdot \mathrm{P}\left[\Delta \ge 0\ge V^{1,i+1}|F_{i,1}=1\right]\\ +& \mathrm{P}\left[{\varphi }^{1,i+1}<\delta |(\Delta \le V^{1,i+1}<0)\wedge F_{i,1}=1\right]& & \cdot \mathrm{P}\left[\Delta \le V^{1,i+1}<0|F_{i,1}=1\right]\\ +& \mathrm{P}\left[{\varphi }^{1,i+1}<\delta |(V^{1,i+1}\le \Delta <0)\wedge F_{i,1}=1\right]& & \cdot \mathrm{P}\left[V^{1,i+1}\le \Delta <0|F_{i,1}=1\right]\\ +& \mathrm{P}\left[{\varphi }^{1,i+1}<\delta |(\Delta \le 0\le V^{1,i+1})\wedge F_{i,1}=1\right]& & \cdot \mathrm{P}\left[\Delta \le 0\le V^{1,i+1}|F_{i,1}=1\right]\end{array}$$

$$\begin{array}{rlrl}=& \mathrm{P}\left[2\cdot V^{1,i+1}-\Delta <\delta |(V^{1,i+1}\ge \Delta \ge 0)\wedge F_{i,1}=1\right]& & \cdot \mathrm{P}\left[V^{1,i+1}\ge \Delta \ge 0|F_{i,1}=1\right]\\ +& \mathrm{P}\left[\Delta <\delta |(\Delta \ge V^{1,i+1}\ge 0)\wedge F_{i,1}=1\right]& & \cdot \mathrm{P}\left[\Delta \ge V^{1,i+1}\ge 0|F_{i,1}=1\right]\\ +& \mathrm{P}\left[-\Delta <\delta |(\Delta \le V^{1,i+1}<0)\wedge F_{i,1}=1\right]& & \cdot \mathrm{P}\left[\Delta \le V^{1,i+1}<0|F_{i,1}=1\right]\\ +& \mathrm{P}\left[\Delta -2\cdot V^{1,i+1}<\delta |(V^{1,i+1}\le \Delta <0)\wedge F_{i,1}=1\right]& & \cdot \mathrm{P}\left[V^{1,i+1}\le \Delta <0|F_{i,1}=1\right]\\ +& \mathrm{P}\left[\Delta -2\cdot V^{1,i+1}<\delta |(\Delta \ge 0\ge V^{1,i+1})\wedge F_{i,1}=1\right]& & \cdot \mathrm{P}\left[\Delta \ge 0\ge V^{1,i+1}|F_{i,1}=1\right]\\ +& \mathrm{P}\left[2\cdot V^{1,i+1}-\Delta <\delta |(\Delta \le 0\le V^{1,i+1})\wedge F_{i,1}=1\right]& & \cdot \mathrm{P}\left[\Delta \le 0\le V^{1,i+1}|F_{i,1}=1\right]\end{array}$$

$$\displaystyle\begin{aligned} =\frac{1}{2}&\cdot\Bigg{(}\frac{1}{2}\cdot\frac{\delta-|\Delta|}{2\cdot\delta}+1\cdot\frac{|\Delta|}{2\cdot\delta}+1\cdot\frac{|\Delta|}{2\cdot\delta}+\frac{1}{2}\cdot\frac{\delta-|\Delta|}{2\cdot\delta}+\frac{\delta-|\Delta|}{2\cdot\delta}\cdot\frac{1}{2}+\frac{\delta-|\Delta|}{2\cdot\delta}\cdot\frac{1}{2}\Bigg{)}&\\ =\frac{1}{2}&&\qed\end{aligned}$$

$$\mathrm{P}[Y=2]=$$

$$=$$

$$\mathrm{P}[Y=k]=$$

$$=$$

$$\cdot \mathrm{P}\left[\underset{v=0}{\overset{k-2}{\bigwedge }}{F}_{i+\lfloor \frac{n+v}{N}\rfloor ,n+v\mathrm{mod}N}=1\right]$$

$$=$$

$${\rm P}\left[\phi^{n,i+1}\geq\delta\;\middle|\;F_{i,n,d}=0\wedge F_{i-1,n,d}=1\right]\geq\frac{1}{2}\cdot\Big{(}1-\frac{1}{2\chi}\Big{)}\enspace.$$

$${\varphi }^{n,i+1}$$

$$=|\chi \cdot V^{n,i}+(c_1\cdot r+c_2\cdot s)\cdot \Delta |+|\Delta -\chi \cdot V^{n,i}-(c_1\cdot r+c_2\cdot s)\cdot \Delta |$$

$$=|\chi \cdot V^{n,i}+(c_1\cdot r+c_2\cdot s)\cdot \Delta |+|-\chi \cdot V^{n,i}-((c_1\cdot r+c_2\cdot s)-1)\cdot \Delta |$$

$$=|\chi \cdot V^{n,i}+(c_1\cdot r+c_2\cdot s)\cdot \Delta |+|\chi \cdot V^{n,i}+((c_1\cdot r+c_2\cdot s)-1)\cdot \Delta |$$

$$\ge |\chi \cdot V^{n,i}+(c_1\cdot r+c_2\cdot s)\cdot \Delta +\chi \cdot V^{n,i}+((c_1\cdot r+c_2\cdot s)-1)\cdot \Delta |$$

$$=|2\cdot \chi \cdot V^{n,i}+(2\cdot (c_1\cdot r+c_2\cdot s)-1)\cdot \Delta |$$

$$\mathrm{P}\left[{\varphi }^{n,i+1}\ge \delta |F_{i,n}=0\wedge F_{i-1,n}=1\right]$$

$$\ge \mathrm{P}\left[|2\chi \cdot V^{n,i}+(2\cdot (c_1\cdot r+c_2\cdot s)-1)\cdot \Delta |\ge \delta |F_{i-1,n}=1\right]$$

$$=\mathrm{P}\left[|2\chi \cdot V^{n,i}+(2\cdot (c_1\cdot r+c_2\cdot s)-1)\cdot \Delta |\ge \delta |\begin{array}{rl}& V^{n,i}\cdot (2\cdot (c_1\cdot r+c_2\cdot s)-1)\cdot \Delta \ge 0\\ & \wedge F_{i-1,n}=1\end{array}\right]$$

$$\cdot \mathrm{P}\left[V^{n,i}\cdot (2\cdot (c_1\cdot r+c_2\cdot s)-1)\cdot \Delta \ge 0|F_{i-1,n}=1\right]$$

$$+\mathrm{P}\left[|2\chi \cdot V^{n,i}+(2\cdot (c_1\cdot r+c_2\cdot s)-1)\cdot \Delta |\ge \delta |\begin{array}{rl}& V^{n,i}\cdot (2\cdot (c_1\cdot r+c_2\cdot s)-1)\cdot (G-x^{n,i})<0\\ & \wedge F_{i-1,n}=1\end{array}\right]$$

$$\cdot \mathrm{P}\left[V^{n,i}\cdot (2\cdot (c_1\cdot r+c_2\cdot s)-1)\cdot \Delta <0|F_{i-1,n}=1\right]$$

$$\ge \mathrm{P}\left[|2\chi \cdot V^{n,i}|\ge \delta |V^{n,i}\cdot (2\cdot (c_1\cdot r+c_2\cdot s)-1)\cdot \Delta \ge 0\wedge F_{i-1,n}=1\right]$$

$$\cdot \mathrm{P}\left[V^{n,i}\cdot (2\cdot (c_1\cdot r+c_2\cdot s)-1)\cdot \Delta \ge 0|F_{i-1,n}=1\right]$$

$$\displaystyle\geq\underbrace{{\rm P}\left[2\chi\cdot|V^{n,i}|\geq\delta\;\middle|\;F_{i-1,n}=1\right]}_{=:\alpha}\cdot\underbrace{{\rm P}\left[V^{n,i}\cdot(2\cdot(c_{1}\cdot r+c_{2}\cdot s)-1)\cdot\Delta\geq 0\right]}_{=:\beta}\geq\Big{(}1-\frac{1}{2\chi}\Big{)}\cdot\frac{1}{2}$$

$$\mathrm{P}\left[V^{n,i}\cdot (2\cdot (c_1\cdot r+c_2\cdot s)-1)\cdot \Delta \ge 0\right]$$

$$=\mathrm{P}\left[V^{n,i}\cdot (2\cdot (c_1\cdot r+c_2\cdot s)-1)\cdot \Delta \ge 0|V^{n,i}\ge 0\right]\cdot \mathrm{P}\left[V^{n,i}\ge 0\right]$$

$$+\mathrm{P}\left[V^{n,i}\cdot (2\cdot (c_1\cdot r+c_2\cdot s)-1)\cdot \Delta \ge 0|V^{n,i}<0\right]\cdot \mathrm{P}\left[V^{n,i}<0\right]$$

$$=\frac{1}{2}\cdot \left(\mathrm{P}\left[(2\cdot (c_1\cdot r+c_2\cdot s)-1)\cdot \Delta \ge 0\right]+\mathrm{P}\left[(2\cdot (c_1\cdot r+c_2\cdot s)-1)\cdot \Delta <0\right]\right)$$

$$=\frac{1}{2}$$

Full text

Declarations of interest

None.

1 Introduction

Background

Particle Swarm Optimization (PSO) is a meta-heuristic for so-called continuous black box optimization problems, which means that the objective function is not explicitly known in form of a, e. g., closed formula. PSO produces good results in a variety of different real world applications. The classical PSO as first introduced by Eberhart and Kennedy [KE95, EK95] in the year 1995 works in (solution candidates improving) iterations and can be very easily implemented and adapted to the users’ applications, which lead to increased attention not only among computer scientists. In order to further improve the performance, many authors present changes to the original, “plain,” or classical PSO scheme (for exact definitions, see Sec. 2) to improve the quality of the returned solution.

A serious problem of PSO that can be sometimes observed is the phenomenon of premature stagnation, i. e., the convergence of the swarm to a non-optimal solution. This phenomenon has been theoretically addressed by Lehre and Witt in [LW13]. To overcome such stagnation, the authors propose Noisy PSO that adds a “noise” term to the velocity at every move. They prove that for the Noisy PSO started on a certain simple $1$-dimensional objective function the first hitting time of the $\delta$-neighborhood of the global optimum is finite. However, as proved in [SW13a], premature stagnation of classical PSO does not occur at all when the search space is $1$-dimensional and the objective function is continuous, i. e., in the $1$-dimensional case, PSO provably finds at least a local optimum, almost surely (in the well defined sense of probability theory). Furthermore, [SW13a] shows a similar result for a slightly modified PSO in the general $D$-dimensional case (for stagnation-related results regarding the unmodified PSO, see [RSW15]). This slightly modified PSO assigns a small random velocity in solely one dimension only if the so-called potential of the swarm – a fundamental, measurable quantity of the swarm – of all particles in this dimension falls below a certain (arbitrary and small) bound $\delta$. In the following, we call such random velocity particle moves forced steps and the PSO variant f-PSO. The f-PSO provably finds a local optimum almost surely [SW13a]. Monitoring the swarm’s potential and increasing it from time to time by the forced steps is the key ingredient to mathematically proving successful convergence to a (local) optimum. In this paper, we will use the forced steps of f-PSO and the potential for overcoming another problem of heuristics, the question when to stop the algorithm.

Problem and new contribution

An important problem that arises in the context of PSO and other iterative optimization methods (see [SW81, p. 23f] and further papers, referenced in the related work section below) is when to stop the iterations and to return the best found admissible solution. Usually, the process is terminated (i) when an upper limit on the number of iterations is reached, or (ii) when an upper limit on the number of evaluations of the objective function is reached, or (iii) when the chance of achieving significant improvements in further iterations is extremely low. The choice of the mentioned two upper limits obviously depends on the concrete objective function which means that the user has to have detailed knowledge of the objective function and to “intervene” by hand. In (iii), it is desired and advantageous that the algorithm adaptively decides when to stop, so external intervention is not necessary anymore. To achieve this for PSO, many criteria were introduced in the literature. For a short overview, see below. A commonly used criterion is the swarm diameter, i. e., the maximum distance between two particles, as a measure for the expansion and thus the movement capability of the swarm. When PSO is extended with forced steps, this and also other criteria of this kind do not work because the expansion of the overall swarm is forcefully kept above a certain value and can no longer converge like it does for the classical PSO algorithm. Just the globally best position found by the swarm converges to optimum.

The goal of this paper is to characterize the behavior of the f-PSO when it is close to a (local) optimum. Our experiments and mathematical investigations show that the number of forced moves does not only increase significantly when the distance to the next local optimum falls below a certain bound, but that additionally the number of forced moves performed close to a local optimum is independent of the objective function. In particular, we prove by potential arguments that the swarm “pulsates” in a cloud around the best solution candidate found so far. Note that a similar behavior has been described recently in [YSG18] for the classical PSO. Therefore, by measuring the frequency of occurrences of forced moves, this frequency can act as a stopping criterion, so the swarm may come to a self-determined halt. All findings are experimentally supported.

Related work

Adaptive stopping criteria for PSO have been investigated by Zielinski et al. [ZPL05] and by Zielinski and Laur [ZL07]. In these papers, a list of upper limit-based and adaptive termination criteria is presented and in experiments applied to ($2$-dimensional) established benchmark functions and a real-world problem, resp. In [KHL*+*07], Kwok et al. introduce a stopping criterion for PSO based on the rate of improvements found by the swarm in a given time interval. Quality of the termination is assured using the non-parametric sign-test enforcing a low false-positive rate. A further stopping approach due to Ong and Fukushimah [OF15] combines PSO with gene matrices.

To name some work beyond PSO, Ribeiro et al. [RRS11] introduce a stopping criterion that stops the GRASP (Greedy Randomized Adaptive Search Procedures) algorithm when the probability of an improvement is below a threshold under an experimentally fitted normal distribution.

Safe et al. [SCPB04] present a study of various aspects associated with the specification of termination conditions for simple genetic algorithms. In [AK00], Aytug and Koehler introduce a stopping criterion that stops a genetic algorithm when the optimal solution is found with a specified confidence and thus no real further progress can be expected.

In a very general, ground breaking investigation, Solis and Wets [SW81] consider random search in general and also address the question of stopping criteria. In this context, Dorea [Dor90] presents two adaptive stopping criteria.

Stopping criteria in the context of multi-objective optimization is presented by Martí et al. [MGBM16].

Organization of paper

Sec. 2 presents the necessary description of the classical and forced PSO (f-PSO) and the definition of the potential of a swarm. In Sec. 3, we experimentally and mathematically analyze the behavior of f-PSO when it comes close to a (local) optimum. In Sec. 4, based on the mathematical analysis in Sec. 3, we present the new termination criteria and show their practicability by experimental evaluations.

2 Definitions

Due to the wide variety of existing PSO variants, we first state the exact “classical” PSO algorithm on which our work is based.

Definition 1** (Classical PSO Algorithm).**

A swarm $\cal S$ of $N$ particles moves through the $D$-dimens-ional search space $\mathbb{R}^{D}$ with $\mathcal{D}=\{1,\ldots,D\}$ being the set of dimensions. Each particle $n\in{\cal S}$ consists of a position $X^{n}\in\mathbb{R}^{D}$, a velocity $V^{n}\in\mathbb{R}^{D}$ and a local attractor $L^{n}\in\mathbb{R}^{D}$, storing the best position particle $n$ has visited so far. Additionally, the particles of the swarm share information via the global attractor $G\in\mathbb{R}^{D}$, describing the best point any particle has visited so far, i. e., as soon as a particle has performed its move111the particles’ moves are executed sequentially, so there is some arbitrary order of the particles., it possibly updates the global attractor immediately.

The actual movement of the swarm is governed by the following movement equations where $\chi$, $c_{1}$, $c_{2}\in\mathbb{R}^{+}$ are some positive constants to be fixed later, and $r$ and $s$ are drawn u. a. r. from $[0,1]^{D}$ every time the equation is applied.

[TABLE]

Here, $\odot$ denotes entrywise multiplication (Hadamard product). The application of the equations on particle $n$ is called the move of $n$. When all particles have executed their moves, the swarm has executed one iteration.

Now we repeat the definition of a swarm’s potential measuring how close it is to convergence, i. e., we describe a measure for its movement. A swarm with high potential should be more likely to reach search points far away from the current global attractor, while the potential of a converging swarm approaches [math]. These considerations lead to the following definition [SW13b] (the original definition in [SW13a] is more complex just for technical reasons):

Definition 2** (Potential).**

Fix a moment of the computation of swarm $\cal S$. For $d\in\mathcal{D}$, the current potential $\Phi_{d}$ of $\cal S$ in dimension $d$ is

[TABLE]

$\Phi=(\Phi_{1},\ldots,\Phi_{D})$* is the total potential of $\cal S$, and $\phi_{d}^{n}$ is the contribution of particle $n$ to the potential of $\cal S$ in dimension $d$.*

Note that the potential of the swarm has an entry in every dimension. The swarm comes to a halt if $\Phi\to(0,\ldots,0)$. So, if the particles come close to $G$, the swarm may stop even if $G$ is a non-optimal point, an incident that is called (premature) stagnation. The single values of $\Phi$ can be actually computed and, hence, be used for decisions on the swarm. Between two different dimensions, the potential difference might be large, and “transferring” potential from one dimension to another is not possible due to the movement equations. On the other hand, along the same dimension the particles influence each other and can transfer potential from one particle to the other. This is the reason why there is no potential of individual particles, but only their contribution to the potential.

To address the phenomenon of stagnation, [SW15] slightly modified the PSO movement equations from Definition 1 as follows by “recharging” potential and, hence, keeping the swarm moving furthermore.

Definition 3** (f-PSO).**

The modified movement of the swarm is governed by the following movement equations where $\chi$, $c_{1}$, $c_{2}$, $\delta\in\mathbb{R}^{+}$ are some positive constants to be fixed later, $r$, $s$ and $t$ are drawn u. a. r. from $[0,1]$ for every move and dimension of a particle.

[TABLE]

If the forced velocity update applies to a particle, we call its move and the corresponding dimension in the move forced. An iteration of the swarm is called forced if during this iteration at least one particle performs a forced move. The whole method is called f-PSO.

If it is necessary to identify the values at the beginning of iteration $i$, we write $X_{d}^{n,i}$, $V_{d}^{n,i}$ etc.

Algorithm 1 provides a formal and detailed overview over f-PSO. The introduction of forced velocity updates guarantees that the swarm (or more precisely, the global attractor $G$) almost surely does not converge to a non-optimal point, but finds a local optimum [SW13a]. In our analysis and for the experiments, we used the common parameter settings $\chi=0.72984$, $c_{1}=1.49617$, $c_{2}=1.49617$ and $N>1$ as suggested and used in [CK02], which are parameter settings that are widely used in the literature.

3 Behavior of the f-PSO algorithm

The idea of the modification of the classical PSO algorithm is to help the swarm overcome (premature) stagnation. The modification is implemented in the calculation of the new velocity during a move of a particle (see Def. 3 and Lines 1 and 1 of Algorithm 1). The new velocity of a particle in the current dimension $d$ is drawn u. a. r. from ${[-\delta,\delta]}$ when for all particles $n^{\prime}\in\mathcal{S}$ their contributions $\phi_{d}^{n^{\prime}}$ to the potential in dimension $d$ is less than $\delta$, which means that the range in which the swarm can optimize is small in dimension $d$.

The main topic of this section is to examine if at such a situation the particle swarm uses from now on only forced moves or, more desirable, the swarm will recover and continue using regular velocity updates. In [SW15], it is shown that indeed the latter is the case, unless the global attractor $G$ is already in the neighborhood of a local optimum. As in this case almost all moves are forced, we will see in Sec. 4 that this can be translated into a criterion to stop the PSO’s execution.

3.1 Experiments

We first introduce the notion of the forcing frequency to quantify the number of applications of the forced velocity update (Line 1 of Algorithm 1).

Definition 4** (Absolute and relative forcing frequency).**

Let $I$ denote a (time) interval of $|I|$ iterations during a run of the f-PSO algorithm.

• (a)

Let $\sigma(I,d)$ denote the number of times forced velocity updates (Line 1) have been executed in interval $I$ for dimension $d$, $d\in\{1,\ldots,D\}$. $\sigma(I,d)$ is called the absolute forcing frequency in dimension $d$ over interval $I$.

• (b)

$\sigma(I)=\sum_{d=1}^{D}\sigma(I,d)$* is the total number of forced velocity updates (Line 1) in interval $I$ counted over all dimensions. $\sigma(I)$ is called the absolute forcing frequency over interval $I$.*

Analogously, the relative forcing frequency is $\sigma(I,d)/|I|$ and $\sigma(I)/|I|$, resp.

To explore the behavior of the particle swarm with respect to the absolute forcing frequency $\sigma(I,d)$, we performed a series of experiments on the well known benchmark functions Schwefel, Rosenbrock, Rastrigin, H. C. Elliptic and Sphere (for a comprehensive overview of these and many other benchmark functions, see [Hel10, Sec 4.2]). Sphere and H. C. Elliptic will be analyzed in detail in Sec. 3.2. The tests were performed with Raß’ HiPPSO [Raß17], a high precision implementation of PSO, in order to rule out the influence of insufficient computer systems’ precision when applying a forced velocity update with $\delta$ close to the precision of a, e. g., long double variable in C++. The tests were performed with $N=3$ particles, $D=30$ dimensions, $\delta=10^{-7}$ and the well-known swarm parameters already mentioned in Sec. 2.

Schwefel function

Fig. 1 presents the measurements obtained when optimizing the Schwefel function with $D=30$. The absolute forcing frequencies per dimension over the intervals $I_{i}=[0\ldots 50\,000\cdot i]$ and the function values $f(G)$ of the current global attractors $G$ at the end of each interval are depicted. One can see that the absolute forcing frequency (gradient in the figure) is relatively small if $f(G)$ is far away from the (unique) optimum value and increases considerably when $f(G)$ approaches the optimum. Additionally, one can see that the gradient of the absolute forcing frequency, i. e., the relative forcing frequency, tends to be constant for all $i\geq i_{0}$ after some value $i_{0}$. This is a showcase and remains true for most of the other tested benchmark functions.

Rosenbrock function

During the tests on the standard benchmark functions, processing the Rosenbrock function showed a slightly different behavior. Metaphorically speaking, there is a small “banana-shaped” valley in this function that leads from a local optimum to the global optimum. If the attractors are in this valley (what they are quite early), the chance to improve in some dimensions is quite small and improvements can easily be voided by setbacks in other dimensions. If this happens, forced velocity updates will be applied in some dimensions while in the other dimensions the particles use the usual, regular velocity updates, which can be seen clearly in Fig. 2. The chance to get moving into the optimum’s direction again is not zero, but quite small which leads to a long phase of near-stagnation as also can be seen in Fig. 2. Eventually the swarm will re-start moving again, but it may take a long time to do so.

Rastrigin function

Another phenomenon may occur on functions like Rastrigin, where there are many local optima. Here, it might happen that the global attractor $G$ is already close to a (good) local optimum, but the local attractor $L^{n}$ of some particle $n$ is close to a worse local optimum. In such a situation, the potential is governed by the distance between $G$ and $L^{n}$, and the condition for executing a forced velocity update is not satisfied, even though a good local near-optimum value $f(G)$ has already been found. But our experiments showed that this happens hardly ever.

3.2 Experiments on the Sphere Function, Phases, and

Theoretical Analysis

We now analyze the swarm’s final behavior with the help of experiments on Sphere (the results for H. C. Elliptic are similar). These experiments show that the final behavior can be classified into two main phases: the approaching phase and the pulsation phase that in turn has three distinct sub-phases: forced phase, lockout phase, recovery phase. These findings hold for all objective functions, as the swarm is in this final phase a kind of a blind searching cloud.

So, when the global attractor reaches the $\delta$-neighborhood of a local optimum, the number of global attractor updates decreases and the distance the global attractor moves becomes small. At that time, the whole swarm will start to approach the global attractor. At some point during this approaching phase, the local attractors will be close to the global attractor and the swarm will loose much of its potential. As described in [SW15], the introduction of forced velocity updates avoids that the actual positions of the particles approach each other closer than a value of about $\delta$ and the swarm enters a kind of pulsation phase: periodically it contracts and expands. However, the global attractor still converges (in the mathematical “infinite time” sense) to the local optimum. Also, the positions of the local attractors may further contract.

It can be observed that when the attractors are all almost at the same position and the particles converge to this position, the number of dimensions that are forced begins to increase. During the pulsation phase almost all particles are forced and the relative forcing frequency $\sigma(I,d)/|I|$ begins to stagnate at a certain value $\partial\sigma_{d}$ (the gradient in the figures). Hence, also $\sigma(I)/|I|$ begins to stagnate at $\partial\sigma$. When this stagnation is reached, the f-PSO will become a periodic pulsating process and will behave almost like a blind search in a box with side-length of order of magnitude $\delta$ around the global attractor. Hence, the objective function $f$ will become irrelevant and the value $\partial\sigma$ is independent of $f$.

3.2.1 Experimental Identification of parameters influencing

$\partial\sigma$

In order to identify the parameters that influence $\partial\sigma$, a series of experiments were performed with the Sphere function. To study a pure pulsation phase, the global and local attractors were initially all set to $\smash{\vec{0}}$, which is the global optimum of this function. With the swarm parameters $\chi$, $c_{1}$ and $c_{2}$ being fixed, the two crucial remaining parameters that might influence $\partial\sigma$ are

• •

the number $D$ of dimensions and

• •

the number $N$ of particles,

whereas the choice of $\delta$ has no influence at all. By changing $\delta$ the range of $\phi_{d}$ causing a forced step is influenced for a dimension $d\in\mathcal{D}$. However, the random velocity assigned by the forced step in this dimension $d$ is changed in the same magnitude. Therefore the change of $\delta$ has no influence on the absolute forcing frequency as these two effects cancel each other.

Fig. 3 and 4 show the effect of changes in the dimension number $D$ and swarm size $N$, resp., on the stagnation value $\partial\sigma$ for fixed $|I|$.

The use of a forced move in a given dimension depends only on the velocities and positions of all particles relative to the global attractor $G$ in this dimension. The dimensions are independent of each other, therefore the stagnation value (recall that $|I|$ is fixed) for each number $D$ of dimensions does not change with an increased number of dimensions as shown in Fig. 3. The same is not true for varying the number of particles, see Fig. 4. If the length $|I|$ of the interval and the number $D$ of dimensions are fixed, increasing the number of particles increases the number of applications of the movement equations by $|I|\cdot D$ for each new particle. Therefore, more dimensions can be forced in a given interval, but there are also more particles that must have a partial potential $\phi$ less than $\delta$.

3.2.2 Mathematical analysis of the sub-phases for

arbitrary objective functions

We can identify three sub-phases during the pulsation phase which are repeated in a fixed order for a given dimension.

• (i)

Starting with a forced phase, where each particle has its move in this dimension forced,

• (ii)

followed by a lockout phase with the probability of a forced move being zero and

• (iii)

finally a recovery phase in which the particles attempt to converge to the global attractor until one move of a particle gets forced again and the cycle repeats itself.

To fully understand why this leads to a less forced frequency $\sigma(I)$ per particle, see Fig. 4, we have to take a closer look at the three sub-phases of the pulsation phase. As described in Def. 3, a dimension $d\in\mathcal{D}$ is forced during the move of a particle when all particles $n^{\prime}$ including the particle itself have contribution $\phi_{d}^{n^{\prime}}$ less than $\delta$ in this dimension $d$. For $i\in\mathbb{N}$, $n\in\mathcal{S}$, and $d\in\mathcal{D}$ let $F_{i,n,d}$ be the $\{0,1\}$-indicator variable being $1$ iff in iteration $i$ and dimension $d$ particle $n$ has the velocity update forced. The global and the local attractors are constant and in particular, it can be assumed that $G_{d}^{i}=G$ for fixed $G$. We get the following observations regarding a forced move in a given dimension.

The first phase we analyze is the lockout phase. Recall that $\phi_{d}^{n,i}=|V_{d}^{n,i}|+|G-X_{d}^{n,i}|$ is the contribution of particle $n$ to the potential of $\cal S$ in dimension $d$ in iteration $i$.

Lemma 1**.**

If $\phi_{d}^{n,i+1}\geq\delta$ for particle $n$, then $\forall n^{\prime}\in\mathcal{S},n^{\prime}>n:F_{i,n^{\prime},d}=0$ and $\forall n^{\prime}\in\mathcal{S},n^{\prime}\leq n:F_{i+1,n^{\prime},d}=0$.

Proof.

With the attractors being constant, the only way for particle $n$ to have contribution $\phi_{d}^{n}$ less than $\delta$ is to make a swarm move that reduces the sum of the velocity and the distance to the global attractor to less than $\delta$. The earliest time this is possible is in the next iteration. ∎

By Lemma 1, we get the length of the lockout phase as $N$, the number of particles.

For the forced phase, we have to look at the probability that for a particle $n$ with a forced move in dimension $d$, $\phi_{d}^{n}<\delta$ still applies in the next iteration.

Lemma 2**.**

*Let $n\in\mathcal{S},d\in\mathcal{D}$. Then for $n<N:{\rm P}[F_{i,n+1,d}=1\;|\;F_{i,n,d}=1]=\frac{1}{2}$ and

${\rm P}[F_{i+1,1,d}=1\;|\;F_{i,N,d}=1]=\frac{1}{2}$*

Proof.

As each dimension is independent and the same calculations are performed in each dimension, we may restrict our analysis to one dimension $d$ and omit the index $d$ in the following proofs.

W. l. o. g., let $n=1$. Let $\Delta=G-X^{1,i}$ the (signed) distance of particle $1$ to the fixed global attractor. Hence, $\phi^{1,i}=|V^{1,i}|+|\Delta|$.

[TABLE]

With this probability we can now analyze the length of the forced phase.

Lemma 3**.**

Given $F_{i,n,d}=1$, define $Y$ as the number of consecutive particle steps, during which dimension $d$ of the respective particle is forced, starting with particle $n$ at iteration $i$. Then for $k\geq 1$ holds ${\rm P}[Y=k]=\frac{1}{2^{k-1}}$.

Proof.

Again, we omit the dimension $d$ from the variables.

By Lemma 1, we know that if for one particle $|V^{n,i+1}|+|G-X^{n,i+1}|\geq\delta$ holds, the next particle cannot be forced and with Lemma 2, ${\rm P}[F_{i,n+1}=1|F_{i,n}=1]=\frac{1}{2}$. By assumption, ${\rm P}[F_{i,n}=1]=1$. By induction on $k$, we have:

[TABLE]

∎

Finally for the recovery phase we cannot give an exact length, as it depends on the positions of the particles in the search space. The following Lemma 4 and Fig. 6 however give us at least an idea of the behavior of the particle during this phase.

Lemma 4**.**

The probability for a particle $n$ to get a partial potential $\phi_{d}^{n}$ greater than $\delta$ when not forced for the first time after a chain of forced moves is

[TABLE]

Proof.

By assumption, $L^{n}=G$. Again, $d$ is omitted from the variables.

Applying the PSO movement equations leads to

[TABLE]

Hence,

[TABLE]

Under the assumption $F_{i-1,n}=1$, we have $V^{n,i}\sim\mathcal{U}(-\delta,\delta)$ and $|V^{n,i}|\sim\mathcal{U}(0,\delta)$ ($\mathcal{U}$ denotes the continuous uniform distribution). Therefore we get $\alpha\geq 1-\frac{1}{2\chi}$ and can derive $\beta=\frac{1}{2}$ by the following equations.

[TABLE]

∎

Note that the lower bound on the probability in Lemma 4 is negative for $\chi<\frac{1}{2}$. The actual probability depends on the position in the search space. With $\chi>\frac{1}{2}$, the probability for a particle $n$ to get partial potential $\phi^{n}>\delta$ is positive independent of the position of the particle. Therefore, the length of the recovery phase will decrease drastically with smaller $\chi$ values.

We may conclude that for a given dimension we get the following behavior in the pulsation phase.

Theorem 1**.**

Let $\cal S$ be a swarm of particles at a (local) optimum of an arbitrary function $f$, i. e., $\forall n\in{\cal S}:L^{n}=G={\rm Opt}(f)$ with ${\rm Opt}(f)$ being a (local) optimum of the function $f$. The swarm repeats three phases independently for each dimension:

• •

the forced phase of length $k\geq 1$ with probability $\frac{1}{2^{k-1}}$,

• •

the lockout phase of length $N$,

• •

and the recovery phase of length $\ell\geq 0$

leading to a stagnation of $\sigma(I)/|I|\to\partial\sigma$ at a (local) optimum depending on $N$ and $D$.

Starting at one forced particle, the chance for the next $k$ particles to be forced has probability $2^{-k+1}$ as stated in Lemma 3 which is independent of the number of particles. At some particle $u$, this sequence of forced particles ends, and in this moment the forced phase ends. With Lemma 1, at least the following $N$ particles will not be forced which constitutes the lockout phase. The next particle that can be forced is $(u+1)\bmod N=v$, and for that all particles $n$ that move between $u$ and $v$ need to keep partial potential $\phi_{d}^{n}<\delta$. All these particles use the regular swarm movement equations from Def. 3, and, with $G$ being constant, they are independent of each other. As shown in Lemma 4, the chance for a particle $n$ to get partial potential $\phi^{n}_{d}>\delta$ can be quite probable depending on the actual swarm parameters. Hence, with increasing $N$ the chance for a particle $n$ to get partial potential $\phi_{d}^{n}>\delta$ increases exponentially in the order of the probability of the complementary event. The exact probability for this event depends on the position of the particle relative to the global attractor. The same is true for the number of iterations a particle needs to fix its partial potential $\phi$, but with normal swarm behavior the particle starts to converge to the global attractor and to decrease the potential. Fig. 6 shows the distance of the particles to the global attractor at the global optimum. We can see that the particles tend to converge to the attractor with high probability thus ending the recovery phase. Therefore, the particles once again start to get forced and the same pattern repeats itself.

4 Termination criteria

In Sec. 3 above, we discussed the behavior of the particle swarm when the global and local attractors are directly at a local optimum. In [SW15], it is shown that with f-PSO the global attractor converges (in the infinite-time sense) to an (at least local) optimum almost surely. However, this convergence can be very slow. Therefore, the particle swarm should be stopped at some time and either be re-started with, when indicated, different parameter settings, or, if the current global attractor does not satisfy the desired quality, to use a different algorithm. As the swarm does not know when it is close to a local optimum and due to the forced steps, the swarm will just converge up to a $\delta$-neighborhood of the global attractor which in turn approached a (local) optimum. With our observations from Sec. 3 this means that the forcing frequency will exceed a fixed value which can be measured.

In the following, we will define a new stopping criterion that stops the algorithm at time when the f-PSO algorithm reaches this frequency. Note that reaching this forcing frequency is just an indication that the swarm is close to an optimum.

4.1 Full stagnation criterion

By Theorem 1, we know that the absolute and relative forcing frequency become fixed if the swarm would be at a (local) optimum, independent of the objective function. So one can use an arbitrary known function with known optimum to experimentally determine this frequency by starting the swarm in the optimum, as we did in Sec. 3.2 for Sphere.

We will use this observation to describe the first simple termination criterion: we just stop the swarm when the measured absolute forcing frequency begins to stagnate and is close to the fixed forcing frequency around a (local) optimum.

Definition 5**.**

Let $\sigma_{\rm stag}(N,D,\mu)$ be the objective function-independent, fixed absolute forcing frequency around an optimum when intervals are considered with interval length $\mu$, and let $\gamma\in\mathbb{N}$. The criterion “$(\gamma,\mu)$-full stop” is to terminate the f-PSO execution iff $\sigma_{\rm stag}(N,D,\mu)-\sigma(I)\leq\gamma$ for some interval $I$ with $|I|=\mu$ .

We tested $(\gamma,\mu)$-full stop in Line 1 of Algorithm 1, i. e., f-PSO, on our benchmark functions with $N=5$ particles, $D=15$ dimensions, $\delta=10^{-7}$ and interval length $\mu=50\,000$. The standard swarm parameters were used. The stagnation frequency $\sigma_{\rm stag}(N,D,\mu)=318\,350$ was measured on the Sphere function with all particles initialized at the global optimum. We used $\gamma=1350$, so we terminated the execution when $\sigma(I)\geq 317\,000$ for intervals $I$ of length $\mu$. If better solutions would be required, one could reduce $\gamma$. We compared the solutions obtained with f-PSO that used $(\gamma,\mu)$-full stop with the results obtained by f-PSO that run for $15\,000\,000$ iterations. Table 1 shows the median (not the average) and standard deviation on 500 runs of the norm of the gradient of the global attractor on these benchmark functions with and without the use of criterion $(\gamma,\mu)$-full stop. The less the gradient, the better the global attractor.

On simple benchmark functions like Sphere criterion $(\gamma,\mu)$-full stop produces sufficiently good output. We see that the time of termination is significantly earlier than with pre-specified number of iterations and that the difference in the gradient at termination is negligible.

4.2 Partial stagnation criterion

On more complex functions like Rosenbrock and Rastrigin, we can identify two scenarios that interfere with the full stop termination criterion.

The first scenario can be observed on the Rosenbrock function. It may occur that some dimensions are already highly optimized, i. e., close to an optimal coordinate, while others are still far away from optimum coordinates, a situation which leads to late termination. With Rosenbrock when the particles and the attractors are in the valley between the local and the global attractor, the chance to improve on the not yet optimized dimensions is voided by the worsening in the already good dimensions. Therefore the swarm will not change the local (and, thus, the global) attractors and the attractors are far away from each other. In this case the optimized dimensions reached the average stagnation frequency $\sigma_{\rm stag}(N,D,\mu)/D$, but the other dimensions will not get forced much or not at all as shown in Fig. 2.

The second scenario becomes visible on the Rastrigin function. This function has many local optima with a global optimum at $\vec{0}$. Unlike to Rosenbrock there is not much difference in the behavior of the particles in the different dimensions. Here, the problem originates from many local optima with only a small difference in function values. This can lead to the global attractor being in a better local optimum than some of the local attractors. Given this situation the particles have a very small probability to actually improve their local attractor. At this point, the algorithm should terminate, but given the different positions of the global and the local attractors these particles cannot decrease their velocity below the necessary bounds to reach the probability of a forced step at the stagnation frequency. The current $\sigma$ will stay around some value, which is different from the stagnation frequency of each dimension and with a fluctuation larger than the fluctuation at the stagnation frequency. However, note that the probability for such an event is quite small and was not encountered in our experiments.

To account for such situations we formulate an extended termination criterion.

Definition 6**.**

Let $\sigma_{\rm stag}(N,D,\mu)$ be as defined in Def. 5, and let $\kappa\in[1,D]$. The criterion “$(\kappa,\gamma,\mu)$-partial stop” is to terminate the f-PSO execution iff $\sigma(I)\geq\kappa\cdot\dfrac{\sigma_{\rm stag}(N,D,\mu)-\gamma}{D}$ for some interval $I$ with $|I|=\mu$ .

Def. 6 introduces the additional parameter $\kappa$. It can be interpreted as how many dimensions are required to reach their stagnation frequency such that f-PSO is stopped. Actually, if $\kappa=D$, we have $(\gamma,\mu)$-full stop. By this parameter, the user can choose the degree of convergence necessary before terminating. On complex objective functions that take a long time to evaluate, a small $\kappa$ may be favorable leading to earlier termination and, hence, saving a lot of running time for applying further subsequent optimization methods to further improve the solution quality. On simple, easily to be evaluated functions such as Sphere, one may choose $\kappa=D$, which results in $(\gamma,\mu)$-full stop, as already mentioned. With $(\kappa,\gamma,\mu)$-partial stop, both interfering scenarios are accounted for. There is no differentiation between dimensions, so it makes no difference that like in the first scenario some dimensions $d$ are forced such they have reached (or are close to) their stagnation value $\partial\sigma_{d}$, while in other dimensions still almost no forcing takes place. Similarly, it might happen, like in the second scenario, that in all dimensions the frequency is less than the stagnation value, but the overall sum has reached the stagnation value. We executed a first series of experiments similar to the case with the criterion $(\gamma,\mu)$-full stop. Again, $N=5$ particles, $D=15$ dimensions, $\delta=10^{-7}$, interval length $\mu=50\,000$, $\sigma_{\rm stag}(N,D,\mu)=318\,350$, and $\gamma=1350$. Tables 2 and 3 show the results of the experiments. We tested $(\kappa,\gamma,\mu)$-partial stop with $\kappa=2$ and $\kappa=8$, and for the fixed iteration limit, we used $15\,000\,000$ iterations. We repeated the experiments 500 times.

We conclude that in most cases a small value of $\kappa$ is sufficient without a significant loss in the quality of the returned solution. However, the Rosenbrock measurement suggests that on complex objective functions, a small value of $\kappa$ might lead to too early termination.

If there are many particles and/or a small number of dimensions, the neighborhood of a local optimum can be reached way sooner. Therefore, the question comes up whether the interval length can be reduced without losing the quality of $(\kappa,\gamma,\mu)$-partial stop. We repeated our experiments for an interval length of $\mu=5\,000$. Here, $N=5$ particles, $D=15$ dimensions, $\delta=10^{-7}$ and adjusted, $\sigma_{\rm stag}(N,D,\mu)=31\,835$, and $\gamma=135$ were applied. The results are presented in Tables 4 and 5.

As presumed for the simple functions, f-PSO with $(\kappa,\gamma,\mu)$-partial stop terminates significantly earlier with only a small decline in the quality of the returned solution. For Rosenbrock, the difference in the quality of the returned solution with different choices of the parameter $\kappa$ becomes more distinct. Furthermore it should be noted that as shown in Fig. 5 the standard deviation increases drastically with shorter intervals. In practice a small tolerance is used when testing if the stagnation value is reached. Given this standard deviation on small intervals, a tolerance that compensates this standard deviation can lead to a too early termination when the stagnation value is not reached and the swarm would recover from the stagnation. This leads to a tradeoff between interval length and the risk of a termination when the swarm has not reached or is close to a local optimum.

5 Conclusions

This paper focused on developing stopping criteria for the f-PSO algorithms. We showed that at an optimum point the swarm behaves like a blind searching algorithm and that therefore the number of forced steps in a time interval, the absolute forcing frequency, is larger than an objection function-independent number. Thus, reaching this frequency it can be presumed that the swarm is close to an optimum solution, and it can be stopped.

Acknowledgements

We would like to thank Alexander Raß for fruitful discussions. Research funded in parts by the School of Engineering of the University of Erlangen-Nuremberg.