E05SBF : NAG Library, Mark 24

Note: before using this routine, please read the Users' Note for your implementation to check the interpretation of bold italicised terms and other implementation-dependent details.

Note: this routine uses optional parameters to define choices in the problem specification and in the details of the algorithm. If you wish to use default settings for all of the optional parameters, you need only read Sections 1 to 9 of this document. If, however, you wish to reset some or all of the settings please refer to Section 10 for a detailed description of the algorithm and to Section 11 for a detailed description of the specification of the optional parameters.

E05SBF is designed to search for the global minimum or maximum of an arbitrary function, subject to general nonlinear constraints, using Particle Swarm Optimization (PSO). Derivatives are not required, although these may be used by an accompanying local minimization routine if desired. E05SBF is essentially identical to E05SAF, with an expert interface and various additional parameters added; otherwise most parameters are identical. In particular, E05SAF does not handle general constraints.

E05SBF has been designed to be particularly effective on SMP systems, by allowing multiple threads to advance subiterations of the algorithm in a highly asynchronous manner. In doing this, the callback routines supplied to E05SBF are called simultaneously on multiple threads, and therefore must themselves be thread-safe. In particular, the arrays IUSER and RUSER provided are classified as OpenMP shared memory, and as such it is imperative that any operations performed on these arrays are done in an appropriate manner. Failure to ensure thread safety of the provided callback functions may result in unpredictable behaviour, including the callback functions returning completely wrong solutions to E05SBF – invalidating any solution returned.

When using an SMP parallel version of this routine, you must indicate that the callback routines are thread-safe by setting the optional argument SMP Callback Thread Safe before calling E05SBF in a multi-threaded environment. See Section 11.2 for more information on this and other SMP options.

Note: the stochastic method used in E05SBF will not produce repeatable answers when run on multiple threads.

Before calling E05SBF, E05ZKF must be called with OPTSTR set to ‘Initialize = e05sbf’. Optional parameters may also be specified by calling E05ZKF before the call to E05SBF.

E05SBF uses a stochastic method based on Particle Swarm Optimization (PSO) to search for the global optimum of a nonlinear function

F

, subject to a set of bound constraints on the variables, and optionally a set of general nonlinear constraints. In the PSO algorithm (see Section 10), a set of particles is generated in the search space, and advances each iteration to (hopefully) better positions using a heuristic velocity based upon inertia, cognitive memory and global memory. The inertia is provided by a decreasingly weighted contribution from a particles current velocity, the cognitive memory refers to the best candidate found by an individual particle and the global memory refers to the best candidate found by all the particles. This allows for a global search of the domain in question.

Further, this may be coupled with a selection of local minimization routines, which may be called during the iterations of the heuristic algorithm, the interior phase, to hasten the discovery of locally optimal points, and after the heuristic phase has completed to attempt to refine the final solution, the exterior phase. Different options may be set for the local optimizer in each phase.

Without loss of generality, the problem is assumed to be stated in the following form:

\underset{x \in R^{ndim}}{minimize} F (x) subject to ℓ \leq (\begin{matrix} x \\ c (x) \end{matrix}) \leq u,

where the objective

F (x)

is a scalar function,

c (x)

is a vector of scalar constraint functions,

x

is a vector in

R^{ndim}

and the vectors

ℓ \leq u

are lower and upper bounds respectively for the

ndim

variables and

ncon

constraints. Both the objective function and the

ncon

constraints may be nonlinear. Continuity of

F

, and the functions

c (x)

, is not essential. For functions which are smooth and primarily unimodal, faster solutions will almost certainly be achieved by using Chapter E04 routines directly.

For functions which are smooth and multi-modal, gradient dependent local minimization routines may be coupled with E05SBF.

For multi-modal functions for which derivatives cannot be provided, particularly functions with a significant level of noise in their evaluation, E05SBF should be used either alone, or coupled with E04CBF.

For heavily constrained problems, E05SBF should either be used alone, or coupled with E04UCF/E04UCA provided the function and the constraints are sufficiently smooth.

The

ndim

lower and upper box bounds on the variable

x

are included to initialize the particle swarm into a finite hypervolume, although their subsequent influence on the algorithm is user determinable (see the option Boundary in Section 11). It is strongly recommended that sensible bounds are provided for all variables and constraints.

E05SBF may also be used to maximize the objective function, or to search for a feasible point satisfying the simple bounds and general constraints (see the option Optimize).

Due to the nature of global optimization, unless a predefined target is provided, there is no definitive way of knowing when to end a computation. As such several stopping heuristics have been implemented into the algorithm. If any of these is achieved, E05SBF will exit with

IFAIL = 1

, and the parameter INFORM will indicate which criteria was reached. See INFORM for more information.

In addition, you may provide your own stopping criteria through MONMOD, OBJFUN and CONFUN.

E05SAF provides a simpler interface, without the inclusion of general nonlinear constraints.

Gill P E, Murray W and Wright M H (1981) Practical Optimization Academic Press

Kennedy J and Eberhart R C (1995) Particle Swarm Optimization Proceedings of the 1995 IEEE International Conference on Neural Networks 1942–1948

Koh B, George A D, Haftka R T and Fregly B J (2006) Parallel Asynchronous Particle Swarm Optimization International Journal for Numerical Methods in Engineering 67(4) 578–595

Vaz A I and Vicente L N (2007) A Particle Swarm Pattern Search Method for Bound Constrained Global Optimization Journal of Global Optimization 39(2) 197–219 Kluwer Academic Publishers

Note: for descriptions of the symbolic variables, see Section 10.

If on entry

IFAIL = 0

or

- 1

, explanatory error messages are output on the current error message unit (as defined by X04AAF).

Note: E05SBF may return useful information for one or more of the following detected errors or warnings.

Errors or warnings detected by the routine:

If

IFAIL = 0

(or

IFAIL = 2

) or

IFAIL = 1

on exit, a criterion will have been reached depending on user selected options. As with all global optimization software, the solution achieved may not be the true global optimum. Various options allow for either greater search diversity or faster convergence to a (local) optimum (See Sections 10 and 11).

Provided the objective function and constraints are sufficiently well behaved, if a local minimizer is used in conjunction with E05SBF, then it is more likely that the final result will at least be in the near vicinity of a local optimum, and due to the global search characteristics of the particle swarm, this solution should be superior to many other local optima.

Caution should be used in accelerating the rate of convergence, as with faster convergence, less of the domain will remain searchable by the swarm, making it increasingly difficult for the algorithm to detect the basins of attraction of superior local optima. Using the options Repulsion Initialize and Repulsion Finalize described in Section 11 will help to overcome this, by causing the swarm to diverge away from the current optimum once no more local improvement is likely.

On successful exit with guaranteed success,

IFAIL = 0

(or

IFAIL = 2

). This may happen if a Target Objective Value is assigned and is reached by the algorithm at a satisfactorily constrained point. It will also occur if a constrained point is found when

Optimize = CONSTRAINTS

is set.

On successful exit without guaranteed success,

IFAIL = 1

is returned. This will happen if another finalization criterion is achieved without the detection of an error.

In both cases, the value of INFORM provides further information as to the cause of the exit.

The memory used by E05SBF is relatively static throughout. Indeed, most of the memory required is used to store the current particle locations, the cognitive particle memories, the particle velocities and the particle weights. As such, E05SBF may be used in problems with high dimension number (

NDIM > 100

) without the concern of computational resource exhaustion, although the probability of successfully locating the global optimum will decrease dramatically with the increase in dimensionality.

Due to the stochastic nature of the algorithm, the result will vary over multiple runs. This is particularly true if parameters and options are chosen to accelerate convergence at the expense of the global search. However, the option

Repeatability = ON

may be set to initialize the internal random number generator using a preset seed, which will result in identical solutions being obtained.

(For SMP users only) The option

Repeatability = ON

will use preset seeds to initialize the random number generator on each thread, however due to the unpredictable nature of parallel communication, this cannot ensure repeatable results when running on multiple threads, even with SMP Thread Overrun set to force synchronization every iteration.

This example uses a particle swarm to find the global minimum of the

2

dimensional Schwefel function:

\underset{x \in R^{2}}{minimize} f = \sum_{j = 1}^{2} x_{j} \sin (\sqrt{|x_{j}|})

subject to the constraints:

\begin{matrix} 3.0 x_{1} - 2.0 x_{2} < 10.0, \\ - 1.0 < x_{1}^{2} - x_{2}^{2} + 3.0 x_{1} x_{2} < 50000.0, \\ -0.9 < \cos ({(x_{1} / 200)}^{2} + (x_{2} / 100)) < 0.9, \\ -500 \leq x_{1} \leq 500, \\ -500 \leq x_{2} \leq 500 . \end{matrix}

The global optimum has an objective value of

f_{\min} = -731.707

, located at

x = (-394.15, -433.48)

. Only the third constraint is active at this point.

The example demonstrates how to initialize and set the options arrays using E05ZKF, how to query options using E05ZLF, and finally how to search for the global optimum using E05SBF. The problem is solved twice, first using E05SBF alone, and secondly by coupling E05SBF with E04UCF/E04UCA as a dedicated local minimizer. In both cases the default option

Repeatability = ON

is used to produce repeatable solutions.

Note: for users of the NAG Library for SMP & Multicore the following example program does not include the setting of the optional parameter SMP Callback Thread Safe, and as such if run on multiple threads it will issue an error message.

None.

The following pseudo-code describes the algorithm used with the repulsion mechanism.

\begin{array}{l} INITIALIZE & for ​ j = 1, ​ n_{p} \\ x_{j} = R \in U (ℓ_{box}, u_{box}) \\ {\hat{x}}_{j} = \{\begin{cases} R \in U (ℓ_{box}, u_{box}) & Start = COLD \\ {\hat{x}}_{j}^{0} & Start = WARM \end{cases} \\ v_{j} = R \in U ({- V}_{\max}, V_{\max}) \\ {\hat{f}}_{j} = \{\begin{cases} F ({\hat{x}}_{j}) & Start = COLD \\ {\hat{f}}_{j}^{0} & Start = WARM \end{cases} \\ {\hat{e}}_{j} = \{\begin{cases} e ({\hat{x}}_{j}) & Start = COLD \\ {\hat{e}}_{j}^{0} & Start = WARM \end{cases} \\ w_{j} = \{\begin{cases} W_{\max} & Weight Initialize = MAXIMUM \\ W_{ini} & Weight Initialize = INITIAL \\ R \in U (W_{\min}, W_{\max}) & Weight Initialize = RANDOMIZED \end{cases} \\ end for \\ \tilde{x} = \frac{1}{2} (ℓ_{box} + u_{box}) \\ \tilde{f} = F (\tilde{x}) \\ \tilde{e} = e (\tilde{x}) \\ I_{c} = I_{s} = 0 \\ SWARM & while (not finalized), \\ I_{c} = I_{c} + 1 \\ for ​ j = 1, n_{p} \\ x_{j} = BOUNDARY (x_{j}, ℓ_{box}, u_{box}) \\ f_{j} = F (x_{j}) \\ e_{j} = e (x_{j}) \\ if ​ (f_{j} / f_{scale} + ϕ (w_{j}) ‖e_{j}‖ < {\hat{f}}_{j} / f_{scale} + ϕ (w_{j}) ‖{\hat{e}}_{j}‖) \\ {\hat{f}}_{j} = f_{j}, ​ {\hat{x}}_{j} = x_{j} \\ if ​ ((‖e_{j}‖ < ‖\tilde{e}‖) ​ or ​ (‖e_{j}‖ \approx ‖\tilde{e}‖ ​ and ​ f_{j} < \tilde{f})) \\ \tilde{f} = f_{j}, ​ \tilde{x} = x_{j} \\ end for \\ if ​ (new (\tilde{f})) \\ LOCMIN (\tilde{x}, \tilde{f}, \tilde{e}, O_{i}), ​ I_{s} = 0 \\ [see note on repulsion below for code insertion] \\ else \\ I_{s} = I_{s} + 1 \\ for ​ j = 1, n_{p} \\ v_{j} = w_{j} v_{j} + C_{s} D_{1} ({\hat{x}}_{j} - x_{j}) + C_{g} D_{2} (\tilde{x} - x_{j}) \\ x_{j} = x_{j} + v_{j} \\ if ​ (‖x_{j} - \tilde{x}‖ < dtol) \\ reset ​ x_{j}, ​ v_{j}, ​ w_{j}; ​ {\hat{x}}_{j} = x_{j} \\ else \\ update ​ (w_{j}) \\ end for \\ if (target achieved or termination criterion satisfied) \\ finalized = true \\ MONMOD (x_{j}) \\ end \\ LOCMIN (\tilde{x}, \tilde{f}, \tilde{e}, O_{e}) \end{array}

The definition of terms used in the above pseudo-code are as follows.

$n_{p}$	the number of particles, NPAR
$ℓ_{box}$	array of NDIM lower box bounds
$u_{box}$	array of NDIM upper box bounds
$x_{j}$	position of particle $j$
${\hat{x}}_{j}$	best position found by particle $j$
$\tilde{x}$	best position found by any particle
$f_{j}$	$F (x_{j})$
${\hat{f}}_{j}$	$F ({\hat{x}}_{j})$ , best value found by particle $j$
$\tilde{f}$	$F (\tilde{x})$ , best value found by any particle
$e_{k} (x)$	$k$ th (scaled) constraint violation at $x$ , evaluated as $\min (c_{k} (x) - l_{NDIM + k}, 0.0) + \max (c_{k} (x) - u_{NDIM + k}, 0.0)$ ; this may be scaled by the maximum $k$ th constraint found thus far
$e (x)$	the array of NCON constraint violations, $e_{k} (x)$ , for $k = 1, 2, \dots, NCON$ , at a point $x$
$e_{j}$	$e (x_{j})$ , the array of constraint violations evaluated at $x_{j}$
${\hat{e}}_{j}$	$e ({\hat{x}}_{j})$ , the array of constraint violations evaluated at ${\hat{x}}_{j}$
$\tilde{e}$	$e (\tilde{x})$ , the array of constraint violations evaluated at $\tilde{x}$
$v_{j}$	velocity of particle $j$
$w_{j}$	weight on $v_{j}$ for velocity update, decreasing according to Weight Decrease
$V_{\max}$	maximum absolute velocity, dependent upon Maximum Variable Velocity
$I_{c}$	swarm iteration counter
$I_{s}$	iterations since $\tilde{x}$ was updated
$f_{scale}$	objective function scaling defined by the options Constraint Scaling, Objective Scaling and Objective Scale.
$D_{1}$ , $D_{2}$	diagonal matrices with random elements in range $(0, 1)$
$C_{s}$	the cognitive advance coefficient which weights velocity towards ${\hat{x}}_{j}$ , adjusted using Advance Cognitive
$C_{g}$	the global advance coefficient which weights velocity towards $\tilde{x}$ , adjusted using Advance Global
$dtol$	the Distance Tolerance for resetting a converged particle
$R \in U (ℓ_{box}, u_{box})$	an array of random numbers whose $i$ -th element is drawn from a uniform distribution in the range $({ℓ_{box}}_{i}, {u_{box}}_{i})$ , for $i = 1, 2, \dots, NDIM$
$O_{i}$	local optimizer interior options
$O_{e}$	local optimizer exterior options
$ϕ (w_{j})$	a function of $w_{j}$ designed to increasingly weight towards minimizing constraint violations as $w_{j}$ decreases
$LOCMIN (x, f, e, O)$	apply local optimizer using the set of options $O$ using the solution $(x, f, e)$ as the starting point, if used (not default)
MONMOD	monitor progress and possibly modify $x_{j}$
BOUNDARY	apply required behaviour for $x_{j}$ outside bounding box, (see Boundary)
new ( $\tilde{f}$ )	true if $\tilde{x}$ , $\tilde{c}$ , $\tilde{f}$ were updated at this iteration

Additionally a repulsion phase can be introduced by changing from the default values of options Repulsion Finalize (

r_{f}

), Repulsion Initialize (

r_{i}

) and Repulsion Particles (

r_{p}

). If the number of static particles is denoted

n_{s}

then the following can be inserted after the new(

\tilde{f}

) check in the pseudo-code above.

\begin{array}{l} else ​ & if ​ & (n_{s} \geq r_{p} ​ and ​ r_{i} \leq I_{s} \leq r_{i} + r_{f}) \\ LOCMIN (\tilde{x}, \tilde{f}, \tilde{e}, O_{i}) \\ use ​ {- C}_{g} ​ instead of ​ C_{g} ​ in velocity updates \\ if ​ & (I_{s} = r_{i} + r_{f}) \\ I_{s} = 0 \end{array}

The algorithm has been parallelized to allow for a high degree of asynchronicity between threads. Each thread is assigned a static number of the NPAR particles requested, and performs a sub-iteration using these particles and a private copy of

\tilde{x}

. The thread only updates this private copy if a superior solution is found.

Once a thread has completed a sub-iteration, it enters a brief critical section where it compares this private

\tilde{x}

to a globally accessible version. If either is superior, the inferior version is updated and the thread continues into a new sub-iteration.

Parallelizing the algorithm in this way allows for individual threads to continue searching even if other threads are completing sub-iterations in inferior times. The optional argument SMP Thread Overrun allows you to force a synchronization across the team of threads once one thread completes sufficiently more sub-iterations than the slowest thread. In particular, this may be used to force synchronization after every sub-iteration if so desired.

This section can be skipped if you wish to use the default values for all optional parameters, otherwise, the following is a list of the optional parameters available and a full description of each optional parameter is provided in Section 11.1.

For each option, we give a summary line, a description of the optional parameter and details of constraints.

The summary line contains:

All options accept the value ‘DEFAULT’ in order to return single options to their default states.

Keywords and character values are case insensitive, however they must be separated by at least one space.

For E05SBF the maximum length of the parameter CVALUE used by E05ZLF is

15

.

The cognitive advance coefficient,

C_{s}

. When larger than the global advance coefficient, this will cause particles to be attracted toward their previous best positions. Setting

r = 0.0

will cause E05SBF to act predominantly as a local optimizer. Setting

r > 2.0

may cause the swarm to diverge, and is generally inadvisable. At least one of the global and cognitive coefficients must be nonzero.

The global advance coefficient,

C_{g}

. When larger than the cognitive coefficient this will encourage convergence toward the best solution yet found. Values

r \in (0, 1)

will inhibit particles overshooting the optimum. Values

r \in [1, 2)

cause particles to fly over the optimum some of the time. Larger values can prohibit convergence. Setting

r = 0.0

will remove any attraction to the current optimum, effectively generating a Monte–Carlo multi-start optimization algorithm. At least one of the global and cognitive coefficients must be nonzero.

Determines the behaviour if particles leave the domain described by the box bounds. This only affects the general PSO algorithm, and will not pass down to any NAG local minimizers chosen.

This option is only effective in those dimensions for which

BL (i) \neq BU (i)

,

i = 1, 2, \dots, NDIM

.

Determines with respect to which norm the constraint residuals should be constructed. These are automatically scaled with respect to NCON as stated. For the set of (scaled) violations

e

, these may be,

Internally, each constraint violation is scaled with respect to the maximum violation yet achieved for that constraint. This option acts as a ceiling for this scale.

Constraint:

r > 1.0

.

Determines whether to scale the constraints and objective function when constructing the penalty function.

The minimum scaled improvement in the constraint violation for a location to be immediately superior to that in memory, regardless of the objective value.

Constraint:

r > 0.0

.

The maximum scaled violation of the constraints for which a sample particle is considered comparable to the current global optimum. Should this not be exceeded, then the current global optimum will be updated if the value of the objective function of the sample particle is superior.

Activates or deactivates the error exit associated with the inability to completely satisfy all constraints,

IFAIL = 4

. It is advisable to deactivate this option if

IFAIL = 0

on entry and the satisfaction of all constraints is not program critical.

Determines whether distances should be scaled by box widths.

This is the distance,

dtol

between particles and the global optimum which must be reached for the particle to be considered converged, i.e., that any subsequent movement of such a particle cannot significantly alter the global optimum. Once achieved the particle is reset into the box bounds to continue searching.

Constraint:

r > 0.0

.

The parameter defines

ε_{r}

, which is intended to be a measure of the accuracy with which the problem function

F (x)

can be computed. If

r < ε

or

r \geq 1

, the default value is used.

The value of

ε_{r}

should reflect the relative precision of

1 + |F (x)|

; i.e.,

ε_{r}

acts as a relative precision when

|F|

is large, and as an absolute precision when

|F|

is small. For example, if

F (x)

is typically of order

1000

and the first six significant digits are known to be correct, an appropriate value for

ε_{r}

would be

10^{- 6}

. In contrast, if

F (x)

is typically of order

10^{- 4}

and the first six significant digits are known to be correct, an appropriate value for

ε_{r}

would be

10^{- 10}

. The choice of

ε_{r}

can be quite complicated for badly scaled problems; see Chapter 8 of Gill et al. (1981) for a discussion of scaling techniques. The default value is appropriate for most simple functions that are computed with full accuracy. However when the accuracy of the computed function values is known to be significantly worse than full precision, the value of

ε_{r}

should be large enough so that no attempt will be made to distinguish between function values that differ by less than the error inherent in the calculation.

Contracts the box boundaries used by a box constrained local minimizer to,

[β_{l}, β_{u}]

, containing the start point

x

, where

\begin{matrix} \partial_{i} = r \times (u_{i} - ℓ_{i}) \\ β_{l}^{i} = \max (ℓ_{i}, x_{i} - \frac{\partial_{i}}{2}) \\ β_{u}^{i} = \min (u_{i}, x_{i} + \frac{\partial_{i}}{2}), i = 1, \dots, NDIM . \end{matrix}

Smaller values of

r

thereby restrict the size of the domain exposed to the local minimizer, possibly reducing the amount of work done by the local minimizer.

Constraint:

0.0 \leq r \leq 1.0

.

The maximum number of iterations or function evaluations the chosen local minimizer will perform inside (outside) the main loop if applicable. For the NAG minimizers these correspond to:

Minimizer	Parameter/option	Default Interior	Default Exterior
E04CBF	MAXCAL	$NDIM + 10$	$2 \times NDIM + 15$
E04DGF/E04DGA	Iteration Limit	$\max (30, 3 \times NDIM)$	$\max (50, 5 \times NDIM)$
E04UCF/E04UCA	Major Iteration Limit	$\max (10, 2 \times NDIM)$	$\max (30, 3 \times NDIM)$

Unless set, these are functions of the parameters passed to E05SBF.

Setting

i = 0

will disable the local minimizer in the corresponding algorithmic region. For example, setting

Local Interior Iterations = 0

and

Local Exterior Iterations = 30

will cause the algorithm to perform no local minimizations inside the main loop of the algorithm, and a local minimization with upto

30

iterations after the main loop has been exited.

Note: currently E04JYF or E04KZF are restricted to using

400 \times NDIM

and

50 \times NDIM

as function evaluation limits respectively. This applies to both local minimizations inside and outside the main loop. They may still be deactivated in either phase by setting

i = 0

, and may subsequently be reactivated in either phase by setting

i > 0

.

Constraint:

i_{1} \geq 0

,

i_{2} \geq 0

.

This is the tolerance provided to a local minimizer in the interior (exterior) of the main loop of the algorithm.

Constraint:

r_{1} > 0.0

,

r_{2} > 0.0

.

Where applicable, the secondary number of iterations the chosen local minimizer will use inside (outside) the main loop. Currently the relevant default values are:

Minimizer	Parameter/option	Default Interior	Default Exterior
E04UCF/E04UCA	Minor Iteration Limit	$\max (10, 2 \times NDIM)$	$\max (30, 3 \times NDIM)$

Constraint:

i_{1} \geq 0

,

i_{2} \geq 0

.

Allows for a choice of Chapter E04 routines to be used as a coupled, dedicated local minimizer.

The maximum number of evaluations of the objective function. When reached this will return

IFAIL = 1

and

INFORM = 6

.

Constraint:

i > 0

.

The maximum number of complete iterations that may be performed. Once exceeded E05SBF will exit with

IFAIL = 1

and

INFORM = 5

.

Unless set, this adapts to the parameters passed to E05SBF.

Constraint:

i \geq 1

.

The maximum number of iterations without any improvement to the current global optimum. If exceeded E05SBF will exit with

IFAIL = 1

and

INFORM = 4

. This exit will be hindered by setting Maximum Iterations Static Particles to larger values.

Constraint:

i \geq 1

.

The minimum number of particles that must have converged to the current optimum before the routine may exit due to Maximum Iterations Static with

IFAIL = 1

and

INFORM = 4

.

Constraint:

i \geq 0

.

The maximum number of particles that may converge to the current optimum. When achieved, E05SBF will exit with

IFAIL = 1

and

INFORM = 3

. This exit will be hindered by setting ‘Repulsion’ options, as these cause the swarm to re-expand.

Constraint:

i > 0

.

The maximum number of particles that may be reset after converging to the current optimum. Once achieved no further particles will be reset, and any particles within Distance Tolerance of the global optimum will continue to evolve as normal.

Constraint:

i > 0

.

Along any dimension

j

, the absolute velocity is bounded above by

|v_{j}| \leq r \times (u_{j} - ℓ_{j}) = V_{\max}

. Very low values will greatly increase convergence time. There is no upper limit, although larger values will allow more particles to be advected out of the box bounds, and values greater than

4.0

may cause significant and potentially unrecoverable swarm divergence.

Constraint:

r > 0.0

.

The initial scale for the objective function. This will remain fixed if

Objective Scaling = USER

is selected.

The method of (re)scaling applied to the objective function when the routine detects a significant difference between the scale and the global and cognitive memory (

\tilde{f}

and

{\hat{f}}_{j}

). This only has an effect when

NCON > 0

and Constraint Scaling is active.

Determines whether to maximize or minimize the objective function, or ignore the objective and search for a constrained point.

Constraint: if

Optimize = CONSTRAINTS

,

NCON > 0

is required.

Allows for the same random number generator seed to be used for every call to E05SBF.

Repeatability = OFF

is recommended in general.

OFF: The internal generation of random numbers will be nonrepeatable.
ON: The same seed will be used.

The number of iterations performed in a repulsive phase before re-contraction. This allows a re-diversified swarm to contract back toward the current optimum, allowing for a finer search of the near optimum space.

Constraint:

i \geq 2

.

The number of iterations without any improvement to the global optimum before the algorithm begins a repulsive phase. This phase allows the particle swarm to re-expand away from the current optimum, allowing more of the domain to be investigated. The repulsive phase is automatically ended if a superior optimum is found.

Constraint:

i \geq 2

.

The number of particles required to have converged to the current optimum before any repulsive phase may be initialized. This will prevent repulsion before a satisfactory search of the near optimum area has been performed, which may happen for large dimensional problems.

Constraint:

i \geq 0

.

Used to affect the initialization of the routine.

COLD: The random number generators and all initialization data will be generated internally. The variables XBEST, FBEST and CBEST need not be set.
WARM: You must supply the initial best location, function and constraint violation values XBEST, FBEST and CBEST. This option is recommended if you already have a data set you wish to improve upon.

The target standard deviation of the particle distances from the current optimum. Once the standard deviation is below this level, E05SBF will exit with

IFAIL = 1

and

INFORM = 2

. This criterion will be penalized by the use of ‘Repulsion’ options, as these cause the swarm to re-expand, increasing the standard deviation of the particle distances from the best point.

In SMP parallel implementations of E05SBF, the standard deviation will be calculated based only on the particles local to the particular thread that checks for finalization. Considerably fewer particles may be used in this calculation than when the algorithm is run in serial. It is therefore recommended that you provide a smaller value of Swarm Standard Deviation when running in parallel than when running in serial.

Constraint:

r \geq 0.0

.

Activate or deactivate the use of a target value as a finalization criterion. If active, then once the supplied target value for the objective function is found (beyond the first iteration if Target Warning is active) E05SBF will exit with

IFAIL = 0

and

INFORM = 1

. Other than checking for feasibility only (

Optimize = CONSTRAINTS

), this is the only finalization criterion that guarantees that the algorithm has been successful. If the target value was achieved at the initialization phase or first iteration and Target Warning is active, E05SBF will exit with

IFAIL = 2

. This option may take any real value

r

, or the character ON/OFF as well as DEFAULT. If this option is queried using E05ZLF, the current value of

r

will be returned in RVALUE, and CVALUE will indicate whether this option is ON or OFF. The behaviour of the option is as follows:

If you have given a target objective value to reach in

objval

(the value of the optional parameter Target Objective Value),

objsfg

sets your desired safeguarded termination tolerance, for when

objval

is close to zero.

Constraint:

objsfg \geq 2 ε

.

The optional tolerance to a user-specified target value.

Constraint:

r \geq 0.0

.

Activates or deactivates the error exit associated with the target value being achieved before entry into the main loop of the algorithm,

IFAIL = 2

.

Adjusts the level of gradient checking performed when gradients are required. Gradient checks are only performed on the first call to the chosen local minimizer if it requires gradients. There is no guarantee that the gradient check will be correct, as the finite differences used in the gradient check are themselves subject to inaccuracies.

Determines how particle weights decrease.

The initial value of any particle's inertial weight,

W_{ini}

, or the minimum possible initial value if initial weights are randomized. When set, this will override

Weight Initialize = RANDOMIZED

or

MAXIMUM

, and as such these must be set afterwards if so desired.

Constraint:

W_{\min} \leq r \leq W_{\max}

.

Determines how the initial weights are distributed.

INITIAL: All weights are initialized at the initial weight, $W_{ini}$ , if set. If Weight Initial has not been set, this will be the maximum weight, $W_{\max}$ .
MAXIMUM: All weights are initialized at the maximum weight, $W_{\max}$ .
RANDOMIZED: Weights are uniformly distributed in $(W_{\min}, W_{\max})$ or $(W_{ini}, W_{\max})$ if Weight Initial has been set.

The maximum particle weight,

W_{\max}

.

Constraint:

1.0 \geq r \geq W_{\min}

(If

W_{ini}

has been set then

1.0 \geq r \geq W_{ini}

.)

The minimum achievable weight of any particle,

W_{\min}

. Once achieved, no further weight reduction is possible.

Constraint:

0.0 \leq r \leq W_{\max}

(If

W_{ini}

has been set then

0.0 \leq r \leq W_{ini}

.)

Determines how particle weights are re-initialized.

The constant

W_{val}

used with

Weight Decrease = INTEREST

.

Constraint:

0.0 \leq r \leq \frac{1}{3}

.

This section details additional options available to users of the NAG Library for SMP & Multicore. In particular it includes the option SMP Callback Thread Safe, which must be set before calling E05SBF with multiple threads.

Declare that the callback routines you provide are or are not thread safe. In particular, this indicates that access to the shared memory arrays IUSER and RUSER from within your provided callbacks is done in a thread-safe manner. If these arrays are just used to pass constant data, then you may assume they are thread safe. If these are also used for workspace, or passing variable data such as random number generator seeds, then you must ensure these are accessed and updated safely. Whilst this can be done using OpenMP critical sections, we suggest their use is minimized to prevent unnecessary bottlenecks, and that instead individual threads have access to independent subsections of the provided arrays where possible.

Determines how many threads will attempt to locally minimize the best found solution after the routine has exited the main loop.

Determines whether the user-supplied function MONMOD is invoked once every sub-iteration each thread performs, or only once by a single thread after all threads have completed at least one sub-iteration.

Determines how many threads support a particle subswarm. This is an extra collection of particles constrained to search only within a hypercube of edge length

10.0 \times Distance Tolerance

of the best point known to an individual thread. This may improve the number of iterations required to find a provided target, particularly if no local minimizer is in use.

If

i \leq 0

, then this will be disabled on all the threads.

If

i \geq OMP_NUM_THREADS

, then all the threads will support a particle subswarm.

This option provides control over the level of asynchronicity present in a simulation. In particular, a barrier synchronization between all threads is performed if any thread completes

i

sub-iterations more than the slowest thread, causing all threads to be exposed to the current best solution. Allowing asynchronous behaviour does however allow individual threads to focus on different global optimum candidates some of the time, which can inhibit convergence to unwanted sub-optima. It also allows for threads to continue searching when other threads are completing sub-iterations at a slower rate.

If

i < 1

, then the algorithm will force a synchronization between threads at the end of each iteration.

INFORM	Meaning
$< 0$	Exit from a user-supplied subroutine.
0	E05SBF has detected an error and terminated.
1	The provided objective target has been achieved. (Target Objective Value).
2	The standard deviation of the location of all the particles is below the set threshold (Swarm Standard Deviation). If the solution returned is not satisfactory, you may try setting a smaller value of Swarm Standard Deviation, or try adjusting the options governing the repulsive phase (Repulsion Initialize, Repulsion Finalize).
3	The total number of particles converged (Maximum Particles Converged) to the current global optimum has reached the set limit. This is the number of particles which have moved to a distance less than Distance Tolerance from the optimum with regard to the $L^{2}$ norm. If the solution is not satisfactory, you may consider lowering the Distance Tolerance. However, this may hinder the global search capability of the algorithm.
4	The maximum number of iterations without improvement (Maximum Iterations Static) has been reached, and the required number of particles (Maximum Iterations Static Particles) have converged to the current optimum. Increasing either of these options will allow the algorithm to continue searching for longer. Alternatively if the solution is not satisfactory, re-starting the application several times with $Repeatability = OFF$ may lead to an improved solution.
5	The maximum number of iterations (Maximum Iterations Completed) has been reached. If the number of iterations since improvement is small, then a better solution may be found by increasing this limit, or by using the option Local Minimizer with corresponding exterior options. Otherwise if the solution is not satisfactory, you may try re-running the application several times with $Repeatability = OFF$ and a lower iteration limit, or adjusting the options governing the repulsive phase (Repulsion Initialize, Repulsion Finalize).
6	The maximum allowed number of function evaluations (Maximum Function Evaluations) has been reached. As with $INFORM = 5$ , increasing this limit if the number of iterations without improvement is small, or decreasing this limit and running the algorithm multiple times with $Repeatability = OFF$ , may provide a superior result.
7	A feasible point has been found. The objective has not been minimized, although it has been evaluated at the final solutions given in XB and XBEST ( $Optimize = CONSTRAINTS$ ).

(i)	The option Target Objective Value has been set and has been reached at a sufficiently constrained point within the search domain.
(ii)	The option $Optimize = CONSTRAINTS$ has been set, and a sufficiently constrained point has been found within the search domain.

(i)	The target was achieved at the first point sampled by the routine. This will be the mean of the lower and upper bounds.
(ii)	The target may have been achieved at a randomly generated sample point. This will always be a possibility provided that the domain under investigation contains a point with a target objective value which is sufficiently constrained.
(iii)	If the Local Minimizer has been set, then a sample point may have been inside the basin of attraction of a satisfactory point. If this occurs repeatedly when the routine is called, it may imply that the objective is largely unimodal, and that it may be more efficient to use the routine selected as the Local Minimizer directly.
(iv)	If $Start = WARM$ has been set, a particle supplied in XBEST with corresponding cognitive memory supplied in FBEST and CBEST satisfies the required target. If $IFAIL = -1$ on entry, the array index corresponding to the first particle to satisfy the target will be returned in the error message. You should verify that the provided values are correct, and if required set a lower Target Objective Value, or lower Constraint Tolerance.

NAG Library Routine Document

E05SBF

+− Contents

1 Purpose

1.1 Information for users of the NAG Library for SMP & Multicore

2 Specification

3 Description

4 References

5 Parameters

6 Error Indicators and Warnings

7 Accuracy

8 Further Comments

9 Example

9.1 Program Text

9.2 Program Data

9.3 Program Results

10 Algorithmic Details

10.1 Details of SMP parallelization

11 Optional Parameters

11.1 Description of the Optional Parameters

11.2 Description of the SMP optional parameters

NAG Library Routine DocumentE05SBF

+− Contents

1 Purpose

1.1 Information for users of the NAG Library for SMP & Multicore

2 Specification

3 Description

4 References

5 Parameters

6 Error Indicators and Warnings

7 Accuracy

8 Further Comments

9 Example

9.1 Program Text

9.2 Program Data

9.3 Program Results

10 Algorithmic Details

10.1 Details of SMP parallelization

11 Optional Parameters

11.1 Description of the Optional Parameters

11.2 Description of the SMP optional parameters

NAG Library Routine Document

E05SBF