NAG FL Interface
e05sbf (nlp_​pso)

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 10 of this document. If, however, you wish to reset some or all of the settings please refer to Section 11 for a detailed description of the algorithm and to Section 12 for a detailed description of the specification of the optional parameters.

1 Purpose

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 arguments added; otherwise most arguments are identical. In particular, e05saf does not handle general constraints.

2 Specification

Fortran Interface
Subroutine e05sbf ( ndim, ncon, npar, xb, fb, cb, bl, bu, xbest, fbest, cbest, objfun, confun, monmod, iopts, opts, iuser, ruser, itt, inform, ifail)
Integer, Intent (In) :: ndim, ncon, npar
Integer, Intent (Inout) :: iopts(*), iuser(*), ifail
Integer, Intent (Out) :: itt(7), inform
Real (Kind=nag_wp), Intent (In) :: bl(ndim+ncon), bu(ndim+ncon)
Real (Kind=nag_wp), Intent (Inout) :: xbest(ndim,npar), fbest(npar), cbest(ncon,npar), opts(*), ruser(*)
Real (Kind=nag_wp), Intent (Out) :: xb(ndim), fb, cb(ncon)
External :: objfun, confun, monmod
C Header Interface
#include <nag.h>
void  e05sbf_ (const Integer *ndim, const Integer *ncon, const Integer *npar, double xb[], double *fb, double cb[], const double bl[], const double bu[], double xbest[], double fbest[], double cbest[],
void (NAG_CALL *objfun)(Integer *mode, const Integer *ndim, const double x[], double *objf, double vecout[], const Integer *nstate, Integer iuser[], double ruser[]),
void (NAG_CALL *confun)(Integer *mode, const Integer *ncon, const Integer *ndim, const Integer *ldcj, const Integer needc[], const double x[], double c[], double cjac[], const Integer *nstate, Integer iuser[], double ruser[]),
void (NAG_CALL *monmod)(const Integer *ndim, const Integer *ncon, const Integer *npar, double x[], const double xb[], const double *fb, const double cb[], const double xbest[], const double fbest[], const double cbest[], const Integer itt[], Integer iuser[], double ruser[], Integer *inform),
Integer iopts[], double opts[], Integer iuser[], double ruser[], Integer itt[], Integer *inform, Integer *ifail)
The routine may be called by the names e05sbf or nagf_glopt_nlp_pso.
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.

3 Description

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 11), 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:
minimize xRndim Fx   subject to   x cx u ,  
where the objective Fx is a scalar function, cx is a vector of scalar constraint functions, x is a vector in Rndim and the vectors 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 cx, 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 12). 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.

4 References

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

5 Arguments

Note: for descriptions of the symbolic variables, see Section 11.
1: ndim Integer Input
On entry: ndim, the number of dimensions.
Constraint: ndim1.
2: ncon Integer Input
On entry: ncon, the number of constraints, not including box constraints.
Constraint: ncon0.
3: npar Integer Input
On entry: npar, the number of particles to be used in the swarm. Assuming all particles remain within constraints, each complete iteration will perform at least npar function evaluations. Otherwise, significantly fewer objective function evaluations may be performed.
Suggested value: npar=10×ndim.
Constraint: npar5×num_threads, where num_threads is the value returned by the OpenMP environment variable OMP_NUM_THREADS, or num_threads is 1 for a serial version of this routine.
4: xbndim Real (Kind=nag_wp) array Output
On exit: the location of the best solution found, x~, in Rndim.
5: fb Real (Kind=nag_wp) Output
On exit: the objective value of the best solution, f~=Fx~.
6: cbncon Real (Kind=nag_wp) array Output
On exit: the constraint violations of the best solution found, e~=ex~. These may have been deemed to be acceptable given the tolerance and scaling of the constraints. See Sections 11 and 12.
7: blndim+ncon Real (Kind=nag_wp) array Input
8: bundim+ncon Real (Kind=nag_wp) array Input
On entry: bl is , the array of lower bounds, bu is u, the array of upper bounds. The first ndim entries in bl and bu must contain the lower and upper simple (box) bounds of the variables respectively. These must be provided to initialize the sample population into a finite hypervolume, although their subsequent influence on the algorithm is user determinable (see the option Boundary in Section 12).
The next ncon entries must contain the lower and upper bounds for any general constraints respectively.
If bli=bui for any i1,,ndim, variable i will remain locked to bli regardless of the Boundary option selected.
It is strongly advised that you place sensible lower and upper bounds on all variables and constraints, even if your model allows for unbounded variables or constraints.
Constraints:
  • blibui, for i=1,2,,ndim+ncon;
  • blibui for at least one i1,,ndim.
9: xbestndimnpar Real (Kind=nag_wp) array Input/Output
Note: the ith component of the best position of the jth particle, x^ji, is stored in xbestij.
On entry: if using Start=WARM, the initial particle positions, x^j0.
On exit: the best positions found, x^j, by the npar particles in the swarm.
10: fbestnpar Real (Kind=nag_wp) array Input/Output
On entry: if using Start=WARM, objective function values, f^j0=Fx^j0, corresponding to the npar particle locations stored in xbest.
On exit: objective function values, f^j=Fx^j, corresponding to the locations returned in xbest.
11: cbestnconnpar Real (Kind=nag_wp) array Input/Output
Note: the kth constraint violation of the jth particle is stored in cbestkj.
On entry: if using Start=WARM, the initial constraint violations, e^j0=ex^j0, corresponding to the npar particle locations.
On exit: the final constraint violations, e^j, corresponding to the locations returned in xbest.
12: objfun Subroutine, supplied by the user. External Procedure
objfun must, depending on the value of mode, calculate the objective function and/or calculate the gradient of the objective function for a ndim-variable vector x. Gradients are only required if a local minimizer has been chosen which requires gradients. See the option Local Minimizer for more information.
The specification of objfun is:
Fortran Interface
Subroutine objfun ( mode, ndim, x, objf, vecout, nstate, iuser, ruser)
Integer, Intent (In) :: ndim, nstate
Integer, Intent (Inout) :: mode, iuser(*)
Real (Kind=nag_wp), Intent (In) :: x(ndim)
Real (Kind=nag_wp), Intent (Inout) :: objf, vecout(ndim), ruser(*)
C Header Interface
void  objfun_ (Integer *mode, const Integer *ndim, const double x[], double *objf, double vecout[], const Integer *nstate, Integer iuser[], double ruser[])
1: mode Integer Input/Output
On entry: indicates which functionality is required.
mode=0
Fx should be returned in objf. The value of objf on entry may be used as an upper bound for the calculation. Any expected value of Fx that is greater than objf may be approximated by this upper bound; that is objf can remain unaltered.
mode=1
Local Minimizer=e04ucf only
First derivatives can be evaluated and returned in vecout. Any unaltered elements of vecout will be approximated using finite differences.
mode=2
Local Minimizer=e04ucf only
Fx must be calculated and returned in objf, and available first derivatives can be evaluated and returned in vecout. Any unaltered elements of vecout will be approximated using finite differences.
mode=5
Fx must be calculated and returned in objf. The value of objf on entry may not be used as an upper bound.
mode=6
Local Minimizer=e04dgf or e04kzf only
All first derivatives must be evaluated and returned in vecout.
mode=7
Local Minimizer=e04dgf or e04kzf only
Fx must be calculated and returned in objf, and all first derivatives must be evaluated and returned in vecout.
On exit: if the value of mode is set to be negative, e05sbf will exit as soon as possible with ifail=3 and inform=mode.
2: ndim Integer Input
On entry: the number of dimensions.
3: xndim Real (Kind=nag_wp) array Input
On entry: x, the point at which the objective function and/or its gradient are to be evaluated.
4: objf Real (Kind=nag_wp) Input/Output
On entry: the value of objf passed to objfun varies with the argument mode.
mode=0
objf is an upper bound for the value of Fx, often equal to the best constraint penalised value of Fx found so far by a given particle if the objective function is strictly positive (see Section 11). Only objective function values less than the value of objf on entry will be used further. As such this upper bound may be used to stop further evaluation when this will only increase the objective function value above the upper bound.
mode=1, 2, 5, 6 or 7
objf is meaningless on entry.
On exit: the value of objf returned varies with the argument mode.
mode=0
objf must be the value of Fx. Only values of Fx strictly less than objf on entry need be accurate.
mode=1 or 6
Need not be set.
mode=2, 5 or 7
Fx must be calculated and returned in objf. The entry value of objf may not be used as an upper bound.
5: vecoutndim Real (Kind=nag_wp) array Input/Output
On entry: if Local Minimizer=e04ucf or , the values of vecout are used internally to indicate whether a finite difference approximation is required. See e04ucf/​e04uca.
On exit: the required values of vecout returned to the calling routine depend on the value of mode.
mode=0 or 5
The value of vecout need not be set.
mode=1 or 2
vecout can contain components of the gradient of the objective function F xi for some i=1,2,ndim, or acceptable approximations. Any unaltered elements of vecout will be approximated using finite differences.
mode=6 or 7
vecout must contain the gradient of the objective function F xi for all i=1,2,ndim. Approximation of the gradient is strongly discouraged, and no finite difference approximations will be performed internally (see e04dgf/​e04dga and e04kzf).
6: nstate Integer Input
On entry: nstate indicates various stages of initialization throughout the routine. This allows for permanent global arguments to be initialized the least number of times. For example, you may initialize a random number generator seed.
nstate=3
SMP users only. objfun is called for the first time in a parallel region on a new thread other than the master thread. You may use this opportunity to set up any thread-dependent information in iuser and ruser.
nstate=2
objfun is called for the very first time. You may save computational time if certain data must be read or calculated only once.
nstate=1
objfun is called for the first time by a NAG local minimization routine. You may save computational time if certain data required for the local minimizer need only be calculated at the initial point of the local minimization.
nstate=0
Used in all other cases.
7: iuser* Integer array User Workspace
8: ruser* Real (Kind=nag_wp) array User Workspace
objfun is called with the arguments iuser and ruser as supplied to e05sbf. You should use the arrays iuser and ruser to supply information to objfun.
objfun must either be a module subprogram USEd by, or declared as EXTERNAL in, the (sub)program from which e05sbf is called. Arguments denoted as Input must not be changed by this procedure.
Note: objfun should not return floating-point NaN (Not a Number) or infinity values, since these are not handled by e05sbf. If your code inadvertently does return any NaNs or infinities, e05sbf is likely to produce unexpected results.
13: confun Subroutine, supplied by the NAG Library or the user. External Procedure
confun must calculate any constraints other than the box constraints. If no constraints are required, confun may be the dummy constraint routine e05szm. (e05szm is included in the NAG Library). For information on how a NAG local minimizer will use confun see the documentation for e04uca.
The specification of confun is:
Fortran Interface
Subroutine confun ( mode, ncon, ndim, ldcj, needc, x, c, cjac, nstate, iuser, ruser)
Integer, Intent (In) :: ncon, ndim, ldcj, needc(ncon), nstate
Integer, Intent (Inout) :: mode, iuser(*)
Real (Kind=nag_wp), Intent (In) :: x(ndim)
Real (Kind=nag_wp), Intent (Inout) :: cjac(ldcj,ndim), ruser(*)
Real (Kind=nag_wp), Intent (Out) :: c(ncon)
C Header Interface
void  confun_ (Integer *mode, const Integer *ncon, const Integer *ndim, const Integer *ldcj, const Integer needc[], const double x[], double c[], double cjac[], const Integer *nstate, Integer iuser[], double ruser[])
1: mode Integer Input/Output
On entry: indicates which values must be assigned during each call of confun. Only the following values need be assigned, for each value of k1,,ncon such that needck>0:
mode=0
the constraint values ckx.
mode=1
rows of the constraint Jacobian, ckxix, for i=1,2,,ndim.
mode=2
the constraint values ckx and the corresponding rows of the constraint Jacobian, ckxix, for i=1,2,,ndim.
On exit: may be set to a negative value if you wish to terminate the solution to the current problem. In this case e05sbf will terminate with ifail=3 and inform=mode as soon as possible.
2: ncon Integer Input
On entry: the number of constraints, not including box bounds.
3: ndim Integer Input
On entry: the number of variables.
4: ldcj Integer Input
On entry: the first dimension of the array cjac as declared in the (sub)program from which e05sbf is called.
5: needcncon Integer array Input
On entry: the indices of the elements of c and/or cjac that must be evaluated by confun. If needck>0, the kth element of c, corresponding to the values of the kth constraint, and/or the available elements of the kth row of cjac, corresponding to the derivatives of the kth constraint, must be evaluated at x (see argument mode).
6: xndim Real (Kind=nag_wp) array Input
On entry: x, the vector of variables at which the constraint functions and/or the available elements of the constraint Jacobian are to be evaluated.
7: cncon Real (Kind=nag_wp) array Output
On exit: if needck>0 and mode=0 or 2, ck must contain the value of ckx. The remaining elements of c, corresponding to the non-positive elements of needc, need not be set.
8: cjacldcjndim Real (Kind=nag_wp) array Input/Output
Note: the derivative of the kth constraint with respect to the ith component, ck xi , is stored in cjacki.
On entry: the elements of cjac are set to special values which enable e05sbf to detect whether they are changed by confun.
On exit: if needck>0 and mode=1 or 2, the elements of cjac corresponding to the kth row of the constraint Jacobian should contain the available elements of the vector ck given by
ck= ck x1 , ck x2 ,, ck xn ,  
where ck xi is the partial derivative of the kth constraint with respect to the ith variable, evaluated at the point x; elements of cjac that remain unaltered will be approximated internally using finite differences. The remaining rows of cjac, corresponding to non-positive elements of needc, need not be set.
It must be emphasized that unassigned elements of cjac are not treated as constant; they are estimated by finite differences, at nontrivial expense. An interval for each element of x is computed automatically at the start of the optimization. The automatic procedure can usually identify constant elements of cjac, which are then computed once only by finite differences.
9: nstate Integer Input
On entry: nstate indicates various stages of initialization throughout the routine. This allows for permanent global arguments to be initialized a minimum number of times. For example, you may initialize a random number generator seed. Note that unless the option Optimize=CONSTRAINTS has been set, objfun will be called before confun.
nstate=3
SMP users only. objfun is called for the first time in a parallel region on a new thread other than the master thread. You may use this opportunity to set up any thread-dependent information in iuser and ruser.
nstate=2
confun is called for the very first time. This argument setting allows you to save computational time if certain data must be read or calculated only once.
nstate=1
confun is called for the first time during a NAG local minimization routine. This argument setting allows you to save computational time if certain data required for the local minimizer need only be calculated at the initial point of the local minimization.
nstate=0
Used in all other cases.
10: iuser* Integer array User Workspace
11: ruser* Real (Kind=nag_wp) array User Workspace
confun is called with the arguments iuser and ruser as supplied to e05sbf. You should use the arrays iuser and ruser to supply information to confun.
confun must either be a module subprogram USEd by, or declared as EXTERNAL in, the (sub)program from which e05sbf is called. Arguments denoted as Input must not be changed by this procedure.
Note: confun should not return floating-point NaN (Not a Number) or infinity values, since these are not handled by e05sbf. If your code inadvertently does return any NaNs or infinities, e05sbf is likely to produce unexpected results.
confun should be tested separately before being used in conjunction with e05sbf.
14: monmod Subroutine, supplied by the NAG Library or the user. External Procedure
A user-specified monitoring and modification function. monmod is called once every complete iteration after a finalization check. It may be used to modify the particle locations that will be evaluated at the next iteration. This permits the incorporation of algorithmic modifications such as including additional advection heuristics and genetic mutations. monmod is only called during the main loop of the algorithm, and as such will be unaware of any further improvement from the final local minimization. If no monitoring and/or modification is required, monmod may be the dummy monitoring routine e05sym. (e05sym is included in the NAG Library).
The specification of monmod is:
Fortran Interface
Subroutine monmod ( ndim, ncon, npar, x, xb, fb, cb, xbest, fbest, cbest, itt, iuser, ruser, inform)
Integer, Intent (In) :: ndim, ncon, npar, itt(7)
Integer, Intent (Inout) :: iuser(*), inform
Real (Kind=nag_wp), Intent (In) :: xb(ndim), fb, cb(ncon), xbest(ndim,npar), fbest(npar), cbest(ncon,npar)
Real (Kind=nag_wp), Intent (Inout) :: x(ndim,npar), ruser(*)
C Header Interface
void  monmod_ (const Integer *ndim, const Integer *ncon, const Integer *npar, double x[], const double xb[], const double *fb, const double cb[], const double xbest[], const double fbest[], const double cbest[], const Integer itt[], Integer iuser[], double ruser[], Integer *inform)
1: ndim Integer Input
On entry: the number of dimensions.
2: ncon Integer Input
On entry: the number of constraints.
3: npar Integer Input
On entry: the number of particles.
4: xndimnpar Real (Kind=nag_wp) array Input/Output
Note: the ith component of the jth particle, xji, is stored in xij.
On entry: the npar particle locations, xj, which will currently be used during the next iteration unless altered in monmod.
On exit: the particle locations to be used during the next iteration.
5: xbndim Real (Kind=nag_wp) array Input
On entry: the location, x~, of the best solution yet found.
6: fb Real (Kind=nag_wp) Input
On entry: the objective value, f~=Fx~, of the best solution yet found.
7: cbncon Real (Kind=nag_wp) array Input
On entry: the constraint violations, e~=ex~, of the best solution yet found.
8: xbestndimnpar Real (Kind=nag_wp) array Input
Note: the ith component of the position of the jth particle's cognitive memory, x^ji, is stored in xbestij.
On entry: the locations currently in the cognitive memory, x^j, for j=1,2,,npar (see Section 11).
9: fbestnpar Real (Kind=nag_wp) array Input
On entry: the objective values currently in the cognitive memory, Fx^j, for j=1,2,,npar.
10: cbestnconnpar Real (Kind=nag_wp) array Input
Note: the kth constraint violation of the jth particle's cognitive memory is stored in cbestkj.
On entry: the constraint violations currently in the cognitive memory, e^=ex^j, for j=1,2,,npar, evaluated at x^j.
11: itt7 Integer array Input
On entry: iteration and function evaluation counters (see description of itt below).
12: iuser* Integer array User Workspace
13: ruser* Real (Kind=nag_wp) array User Workspace
monmod is called with the arguments iuser and ruser as supplied to e05sbf. You should use the arrays iuser and ruser to supply information to monmod.
14: inform Integer Input/Output
On entry: inform=thread_num, where thread_num is the value returned by a call of x06adf. If running in serial this will always be zero.
On exit: setting inform<0 will cause near immediate exit from e05sbf. This value will be returned as inform with ifail=3. You need not set inform unless you wish to force an exit.
monmod must either be a module subprogram USEd by, or declared as EXTERNAL in, the (sub)program from which e05sbf is called. Arguments denoted as Input must not be changed by this procedure.
15: iopts* Integer array Communication Array
Note: the dimension of this array is dictated by the requirements of associated functions that must have been previously called. This array must be the same array passed as argument iopts in the previous call to e05zkf.
On entry: optional parameter array as generated and possibly modified by calls to e05zkf. The contents of iopts must not be modified directly between calls to e05sbf, e05zkf or e05zlf.
16: opts* Real (Kind=nag_wp) array Communication Array
Note: the dimension of this array is dictated by the requirements of associated functions that must have been previously called. This array must be the same array passed as argument opts in the previous call to e05zkf.
On entry: optional parameter array as generated and possibly modified by calls to e05zkf. The contents of opts must not be modified directly between calls to e05sbf, e05zkf or e05zlf.
17: iuser* Integer array User Workspace
iuser is not used by e05sbf, but is passed directly to objfun, confun and monmod and may be used to pass information to these routines.
With care, you may also write information back into iuser. This might be useful, for example, should there be a need to preserve the state of a random number generator.
With SMP-enabled versions of e05sbf the array iuser provided are classified as OpenMP shared memory. Use of iuser has to take account of this in order to preserve thread safety whenever information is written back to either of these arrays.
18: ruser* Real (Kind=nag_wp) array User Workspace
ruser is not used by e05sbf, but is passed directly to objfun, confun and monmod and may be used to pass information to these routines.
With care, you may also write information back into ruser. This might be useful, for example, should there be a need to preserve the state of a random number generator.
With SMP-enabled versions of e05sbf the array ruser provided are classified as OpenMP shared memory. Use of ruser has to take account of this in order to preserve thread safety whenever information is written back to either of these arrays.
19: itt7 Integer array Output
On exit: integer iteration counters for e05sbf.
itt1
Number of complete iterations.
itt2
Number of complete iterations without improvement to the current optimum.
itt3
Number of particles converged to the current optimum.
itt4
Number of improvements to the optimum.
itt5
Number of function evaluations performed.
itt6
Number of particles reset.
itt7
Number of violated constraints at completion. Note this is always calculated using the L1 norm and a nonzero result does not necessarily mean that the algorithm did not find a suitably constrained point with respect to the single norm used.
20: inform Integer Output
On exit: indicates which finalization criterion was reached. The possible values of inform are:
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 L2 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).
If you wish to continue from the final position gained from a previous simulation with adjusted options, you may set the option Start=WARM, and pass back in the returned arrays xbest, fbest, and cbest. You should either record the returned values of xb, fb and cb for comparison, as these will not be re-used by the algorithm, or include them in xbest, fbest and cbest respectively by overwriting the entries corresponding to one particle with the relevant information.
21: ifail Integer Input/Output
On entry: ifail must be set to 0, -1 or 1. If you are unfamiliar with this argument you should refer to Section 7 in the Introduction to the NAG Library CL Interface for details.
On exit: the most common exit will be ifail=1.
For this reason, the value -1 or 1 is recommended. If the output of error messages is undesirable, the value 1 is recommended; otherwise, the recommended value is -1. When the value -1 or 1 is used it is essential to test the value of ifail on exit.
e05sbf returns ifail=0 if and only if a finalization criterion has been reached which can guarantee success. This may only happen if:
  1. (i)The option Target Objective Value has been set and has been reached at a sufficiently constrained point within the search domain.
  2. (ii)The option Optimize=CONSTRAINTS has been set, and a sufficiently constrained point has been found within the search domain.
These finalization criteria are not active using default option settings, and must be explicitly set using e05zkf if required.
e05sbf will return ifail=1 if no error has been detected, and a finalization criterion has been achieved which cannot guarantee success. This does not indicate that the routine has failed, merely that the returned solution cannot be guaranteed to be the true global optimum.
The value of inform should be examined to determine which finalization criterion was reached.
Other positive values of ifail indicate that either an error or a warning has been triggered. See Sections 6, 7 and 11 for more information.

6 Error Indicators and Warnings

If on entry ifail=0 or -1, explanatory error messages are output on the current error message unit (as defined by x04aaf).
Errors or warnings detected by the routine:
ifail=1
A finalization criterion was reached that cannot guarantee success.
On exit, inform=value.
ifail=2
If the option Target Warning has been activated, this indicates that the Target Objective Value has been achieved to specified tolerances at a sufficiently constrained point, either during the initialization phase, or during the first two iterations of the algorithm. While this is not necessarily an error, it may occur if:
  1. (i)The target was achieved at the first point sampled by the routine. This will be the mean of the lower and upper bounds.
  2. (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.
  3. (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.
Assuming that objfun is correct, you may wish to set a better Target Objective Value, or a stricter Target Objective Tolerance.
ifail=3
User requested exit value during call to confun.
User requested exit value during call to monmod.
User requested exit value during call to objfun.
ifail=4
Unable to locate strictly feasible point. value constraints remain violated. This exit may be suppressed using the option Constraint Warning.
ifail=11
On entry, ndim=value.
Constraint: ndim1.
ifail=12
On entry, npar=value.
Constraint: npar5×num_threads, where num_threads is the value returned by the OpenMP environment variable OMP_NUM_THREADS, or num_threads is 1 for a serial version of this routine.
ifail=13
On entry, ncon=value.
Constraint: ncon0.
ifail=14
On entry, blvalue=value and buvalue=value.
Constraint: buibli for all i.
On entry, bli=bui for all box bounds i.
Constraint: bui>bli for at least one box bound i.
ifail=17
e05sbf has been called with ncon>0 and the dummy constraint function e05szm. Only use e05szm with ncon=0.
ifail=18
The option Optimize=CONSTRAINTS is active, however ncon=0.
ifail=19
Error value occurred whilst adjusting to exterior local minimizer options.
Error value occurred whilst adjusting to interior local minimizer options.
ifail=21
Either the option arrays have not been initialized for e05sbf, or they have become corrupted.
ifail=32
Derivative checks indicate possible errors in the supplied derivatives. Gradient checks may be disabled by setting Verify Gradients=OFF.
ifail=51
Multiple SMP threads have been detected; however, the option SMP Callback Thread Safe has not been set.
Set SMP Callback Thread Safe=YES if the provided callbacks are thread safe.
Set SMP Callback Thread Safe=NO if the provided callbacks are not thread safe, to force serial execution.
ifail=-99
An unexpected error has been triggered by this routine. Please contact NAG.
See Section 7 in the Introduction to the NAG Library FL Interface for further information.
ifail=-399
Your licence key may have expired or may not have been installed correctly.
See Section 8 in the Introduction to the NAG Library FL Interface for further information.
ifail=-999
Dynamic memory allocation failed.
See Section 9 in the Introduction to the NAG Library FL Interface for further information.

7 Accuracy

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 11 and 12).
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 12 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.

8 Parallelism and Performance

e05sbf is threaded by NAG for parallel execution in multithreaded implementations of the NAG Library.
e05sbf makes calls to BLAS and/or LAPACK routines, which may be threaded within the vendor library used by this implementation. Consult the documentation for the vendor library for further information.
Please consult the X06 Chapter Introduction for information on how to control and interrogate the OpenMP environment used within this routine. Please also consult the Users' Note for your implementation for any additional implementation-specific information.
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 x~. The thread only updates this private copy if a superior solution is found. In these implementations, this routine may make calls to the user-supplied functions from within an OpenMP parallel region. Thus OpenMP directives within the user functions can only be used if you are compiling the user-supplied function and linking the executable in accordance with the instructions in the Users' Note for your implementation.
Once a thread has completed a sub-iteration, it enters a brief critical section where it compares this private 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 parameter 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.
When using an SMP parallel version of this routine, you must indicate that the callback routines are thread safe by setting the optional parameter SMP Callback Thread Safe before calling e05sbf in a multithreaded environment. See Section 12.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.

9 Further Comments

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 arguments 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.

10 Example

This example uses a particle swarm to find the global minimum of the two-dimensional Schwefel function:
minimize xR2 f = j=1 2 xj sinxj  
subject to the constraints:
3.0 x1 - 2.0 x2 <10.0 , -1.0 < x12 - x22 + 3.0 x1 x2 < 50000.0 , -0.9 < cos x1 / 200 2 + x2 / 100 < 0.9 , -500 x1 500 , -500 x2 500 .  
The global optimum has an objective value of fmin=-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 multithreaded implementations of the NAG Library 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. See the additional example program provided for e05saf for more information on how to safely access independent subsections of the provided iuser and ruser arrays from multiple threads and how to use e05zkf to set additional SMP threading related options.

10.1 Program Text

Program Text (e05sbfe.f90)

10.2 Program Data

None.

10.3 Program Results

Program Results (e05sbfe.r)

11 Algorithmic Details

The following pseudo-code describes the algorithm used with the repulsion mechanism.
INITIALIZE for ​ j=1, ​ np xj = R Ubox,ubox x^ j = R Ubox,ubox Start=COLD x^ j 0 Start=WARM vj = R U -V max ,Vmax f^j = F x^ j Start=COLD f^ j 0 Start=WARM e^ j = e x^ j Start=COLD e^ j 0 Start=WARM wj = Wmax Weight Initialize=MAXIMUM Wini Weight Initialize=INITIAL R U W min , W max Weight Initialize=RANDOMIZED end for x~ = 12 box + ubox f~ = F x~ e~ = e x~ Ic = Is = 0 SWARM while (not finalized), Ic = Ic + 1 for ​ j = 1 , np xj = BOUNDARYxj,box,ubox fj = F xj ej = e xj if ​ fj / fscale + ϕ wj ej < f^j / fscale + ϕwj e^j f^j = fj , ​ x^ j = xj if ​ e j < e~ ​ or ​ e j e~ ​ and ​ fj < f~ f~ = fj , ​ x~ = xj end for if ​ new f~ LOCMINx~,f~,e~,Oi , ​ Is=0 [see note on repulsion below for code insertion] else Is = Is + 1 for ​ j = 1 , np vj = wj vj + Cs D1 x^j - xj + Cg D2 x~ - xj xj = xj + vj if ​ xj - x~ < dtol reset ​ xj, ​ vj, ​ wj; ​ x^j = xj else update ​wj end for if (target achieved or termination criterion satisfied) finalized=true monmod xj end LOCMINx~,f~,e~,Oe  
The definition of terms used in the above pseudo-code are as follows.
np the number of particles, npar
box array of ndim lower box bounds
ubox array of ndim upper box bounds
xj position of particle j
x^j best position found by particle j
x~ best position found by any particle
fj Fxj
f^j Fx^j, best value found by particle j
f~ Fx~, best value found by any particle
ekx kth (scaled) constraint violation at x, evaluated as
minckx-lndim+k,0.0+maxckx-undim+k,0.0; this may be scaled by the maximum kth constraint found thus far
ex the array of ncon constraint violations, ekx, for k=1,2,,ncon, at a point x
ej exj, the array of constraint violations evaluated at xj
e^j ex^j, the array of constraint violations evaluated at x^j
e~ ex~, the array of constraint violations evaluated at x~
vj velocity of particle j
wj weight on vj for velocity update, decreasing according to Weight Decrease
Vmax maximum absolute velocity, dependent upon Maximum Variable Velocity
Ic swarm iteration counter
Is iterations since x~ was updated
fscale objective function scaling defined by the options Constraint Scaling, Objective Scaling and Objective Scale.
D1,D2 diagonal matrices with random elements in range 0,1
Cs the cognitive advance coefficient which weights velocity towards x^j, adjusted using Advance Cognitive
Cg the global advance coefficient which weights velocity towards x~, adjusted using Advance Global
dtol the Distance Tolerance for resetting a converged particle
RUbox,ubox an array of random numbers whose ith element is drawn from a uniform distribution in the range boxi,uboxi, for i=1,2,,ndim
Oi local optimizer interior options
Oe local optimizer exterior options
ϕwj a function of wj designed to increasingly weight towards minimizing constraint violations as wj decreases
LOCMINx,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 xj
BOUNDARY apply required behaviour for xj outside bounding box, (see Boundary)
new (f~) true if x~, c~, f~ were updated at this iteration
Additionally a repulsion phase can be introduced by changing from the default values of options Repulsion Finalize (rf), Repulsion Initialize (ri) and Repulsion Particles (rp). If the number of static particles is denoted ns then the following can be inserted after the new(f~) check in the pseudo-code above.
else ​ if ​ ( ns rp ​ and ​ ri Is ri + rf ) LOCMINx~,f~,e~,Oi use ​ -Cg ​ instead of ​ Cg ​ in velocity updates if ​ Is = ri + rf Is = 0  

12 Optional Parameters

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 12.1.

12.1 Description of the Optional Parameters

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 argument cvalue used by e05zlf is 15.
Advance CognitiverDefault =2.0
The cognitive advance coefficient, Cs. 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.
Advance GlobalrDefault =2.0
The global advance coefficient, Cg. When larger than the cognitive coefficient this will encourage convergence toward the best solution yet found. Values r0,1 will inhibit particles overshooting the optimum. Values r1,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.
BoundaryaDefault =FLOATING
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 blibui, i=1,2,,ndim.
IGNORE
The box bounds are ignored. The objective function is still evaluated at the new particle position.
RESET
The particle is re-initialized inside the domain. x^j, f^j and e^j are not affected.
FLOATING
The particle position remains the same, however the objective function will not be evaluated at the next iteration. The particle will probably be advected back into the domain at the next advance due to attraction by the cognitive and global memory.
HYPERSPHERICAL
The box bounds are wrapped around an ndim-dimensional hypersphere. As such a particle leaving through a lower bound will immediately re-enter through the corresponding upper bound and vice versa. The standard distance between particles is also modified accordingly.
FIXED
The particle rests on the boundary, with the corresponding dimensional velocity set to 0.0.
Constraint NormaDefault =L1
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,
L1
The L1 norm will be used, e1=1ncon1nconek
L2
The L2 norm will be used, e2=1ncon1nconek2
L2SQ
The square of the L2 norm will be used, e22=1ncon1nconek2
LMAX
The L norm will be used, e=max0<knconek
Constraint Scale MaximumrDefault =1.0E6
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.
Constraint ScalingaDefault =INITIAL
Determines whether to scale the constraints and objective function when constructing the penalty function.
OFF
Neither the constraint violations nor the objective will be scaled automatically. This should only be used if the constraints and objective are similarly scaled everywhere throughout the domain.
INITIAL
The maximum of the initial cognitive memories, f^j and e^j, will be used to scale the objective function and constraint violations respectively.
ADAPTIVE
Initially, the maximum of the initial cognitive memories, f^j and e^j, will be used to scale the objective function and constraint violations respectively. If a significant change is detected in the behaviour of the constraints or the objective, these will be rescaled with respect to the current state of the cognitive memory.
Constraint SuperiorityrDefault =0.01
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.
Constraint TolerancerDefault =10-4
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.
Constraint WarningaDefault =ON
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.
OFF
No error will be returned.
ON
An error will be returned if any constraints are sufficiently violated at the end of the simulation.
Distance ScalingaDefault =ON
Determines whether distances should be scaled by box widths.
ON
When a distance is calculated between x and y, a scaled L2 norm is used.
L2 x,y = i | ui i , indim xi - yi ui - i 2 1 2 .  
OFF
Distances are calculated as the standard L2 norm without any rescaling.
L2 x,y = i=1 ndim xi - yi 2 1 2 .  
Distance TolerancerDefault =10-4
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.
Function PrecisionrDefault =ε0.9
The parameter defines εr, which is intended to be a measure of the accuracy with which the problem function Fx can be computed. If r<ε or r1, the default value is used.
The value of εr should reflect the relative precision of 1+Fx; i.e., εr acts as a relative precision when F is large, and as an absolute precision when F is small. For example, if Fx 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 Fx 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.
Local Boundary RestrictionrDefault =0.5
Contracts the box boundaries used by a box constrained local minimizer to, βl,βu, containing the start point x, where
i = r × ui - i βli = maxi, xi - i2 βui = minui, xi + i2 ,  i=1,,ndim .  
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.0r1.0.
Local Interior Iterationsi1
Local Interior Major Iterationsi1
Local Exterior Iterationsi2
Local Exterior Major Iterationsi2
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×ndim+15
e04dgf/​e04dga Iteration Limit max30,3×ndim max50,5×ndim
e04ucf/​e04uca Major Iteration Limit max10,2×ndim max30,3×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×ndim and 50×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: i10i20.
Local Interior Tolerancer1Default =10-4
Local Exterior Tolerancer2Default =10-4
This is the tolerance provided to a local minimizer in the interior (exterior) of the main loop of the algorithm.
Constraint: r1>0.0,r2>0.0.
Local Interior Minor Iterationsi1
Local Exterior Minor Iterationsi2
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 max10,2×ndim max30,3×ndim
Constraint: i10i20.
Local MinimizeraDefault =OFF
Allows for a choice of Chapter E04 routines to be used as a coupled, dedicated local minimizer.
Local Minimizer=OFF
No local minimization will be performed in either the INTERIOR or EXTERIOR sections of the algorithm.
Local Minimizer=e04cbf
Use e04cbf as the local minimizer. This does not require the calculation of derivatives.
On a call to objfun during a local minimization, mode=5.
Local Minimizer=e04kzf
Use e04kzf as the local minimizer. This requires the calculation of derivatives in objfun, as indicated by mode.
The box bounds forwarded to this routine from e05sbf will have been acted upon by Local Boundary Restriction. As such, the domain exposed may be greatly smaller than that provided to e05sbf.
Accurate derivatives must be provided to this routine, and will not be approximated internally. Each iteration of this local minimizer also requires the calculation of both the objective function and its derivative. Hence on a call to objfun during a local minimization, mode=7.
Local Minimizer=e04jyf
Use e04jyf as the local minimizer. This does not require the calculation of derivatives.
On a call to objfun during a local minimization, mode=5.
The box bounds forwarded to this routine from e05sbf will have been acted upon by Local Boundary Restriction. As such, the domain exposed may be greatly smaller than that provided to e05sbf.
Local Minimizer=e04dgf
Use e04dgf/​e04dga as the local minimizer.
Accurate derivatives must be provided, and will not be approximated internally. Additionally, each call to objfun during a local minimization will require either the objective to be evaluated alone, or both the objective and its gradient to be evaluated. Hence on a call to objfun, mode=5 or 7.
Local Minimizer=e04ucf
Use e04ucf/​e04uca as the local minimizer. This operates such that any derivatives of either the objective function or the constraint Jacobian, which you cannot supply, will be approximated internally using finite differences.
Either, the objective, objective gradient, or both may be requested during a local minimization, and as such on a call to objfun, mode=1, 2 or 5.
The box bounds forwarded to this routine from e05sbf will have been acted upon by Local Boundary Restriction. As such, the domain exposed may be greatly smaller than that provided to e05sbf.
Maximum Function EvaluationsiDefault =Imax
The maximum number of evaluations of the objective function. When reached this will return ifail=1 and inform=6.
Constraint: i>0.
Maximum Iterations CompletediDefault =1000×ndim
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: i1.
Maximum Iterations StaticiDefault =100
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: i1.
Maximum Iterations Static ParticlesiDefault =0
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: i0.
Maximum Particles ConvergediDefault =Imax
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.
Maximum Particles ResetiDefault =Imax
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.
Maximum Variable VelocityrDefault =0.25
Along any dimension j, the absolute velocity is bounded above by vjr×uj-j=Vmax. 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.
Objective ScalerDefault =1.0
The initial scale for the objective function. This will remain fixed if Objective Scaling=USER is selected.
Objective ScalingaDefault =MAXIMUM
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 (f~ and f^j). This only has an effect when ncon>0 and Constraint Scaling is active.
MAXIMUM
The objective is rescaled with respect to the maximum absolute value of the objective in the cognitive and global memory.
MEAN
The objective is rescaled with respect to the mean absolute value of the objective in the cognitive and global memory.
USER
The scale remains fixed at the value set using Objective Scale.
OptimizeaDefault =MINIMIZE
Determines whether to maximize or minimize the objective function, or ignore the objective and search for a constrained point.
MINIMIZE
The objective function will be minimized.
MAXIMIZE
The objective function will be maximized. This is accomplished by minimizing the negative of the objective.
CONSTRAINTS
The objective function will be ignored, and the algorithm will attempt to find a feasible point given the provided constraints. The objective function will be evaluated at the best point found with regards to constraint violations, and the final positions returned in xbest. The objective will be calculated at the best point found in terms of constraints only. Should a constrained point be found, e05sbf will exit with ifail=0 and inform=6.
Constraint: if Optimize=CONSTRAINTS, ncon>0 is required.
RepeatabilityaDefault =OFF
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.
Repulsion FinalizeiDefault =Imax
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: i2.
Repulsion InitializeiDefault =Imax
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: i2.
Repulsion ParticlesiDefault =0
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: i0.
SeediDefault =0
Sets the random number generator seed to be used when Repeatability=ON. If set to 0, the default seed will be used. If not, the absolute value of Seed will be used to generate the random number generator seed.
StartaDefault =COLD
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.
Swarm Standard DeviationrDefault =0.1
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: r0.0.
Target ObjectiveaDefault =OFF
Target Objective ValuerDefault =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:
r
Once a point is found with an objective value within the Target Objective Tolerance of r, e05sbf will exit successfully with ifail=0 and inform=1.
OFF
The current value of r will remain stored, however it will not be used as a finalization criterion.
ON
The current value of r stored will be used as a finalization criterion.
DEFAULT
The stored value of r will be reset to its default value (0.0), and this finalization criterion will be deactivated.
Target Objective SafeguardrDefault =10.0ε
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: objsfg2ε.
Target Objective TolerancerDefault =0.0
The optional tolerance to a user-specified target value.
Constraint: r0.0.
Target WarningaDefault =OFF
Activates or deactivates the error exit associated with the target value being achieved before entry into the main loop of the algorithm, ifail=2.
OFF
No error will be returned, and the routine will exit normally.
ON
An error will be returned if the target objective is reached prematurely, and the routine will exit with ifail=2.
Verify GradientsaDefault =ON
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.
OFF
No gradient checking will be performed.
ON
A cheap gradient check will be performed on both the gradients corresponding to the objective through objfun and those provided via the constraint Jacobian through confun.
OBJECTIVE
A more expensive gradient check will be performed on the gradients corresponding to the objective objfun. The gradients of the constraints will not be checked.
CONSTRAINTS
A more expensive check will be performed on the elements of cjac provided via confun. The objective gradient will not be checked.
FULL
A more expensive check will be performed on both the gradient of the objective and the constraint Jacobian.
Weight DecreaseaDefault =INTEREST
Determines how particle weights decrease.
OFF
Weights do not decrease.
INTEREST
Weights decrease through compound interest as wIT+1=wIT1-Wval, where Wval is the Weight Value and IT is the current number of iterations.
LINEAR
Weights decrease linearly following wIT+1=wIT-IT×Wmax-Wmin/ITmax, where IT is the iteration number and ITmax is the maximum number of iterations as set by Maximum Iterations Completed.
Weight InitialrDefault =Wmax
The initial value of any particle's inertial weight, Wini, 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: WminrWmax.
Weight InitializeaDefault =MAXIMUM
Determines how the initial weights are distributed.
INITIAL
All weights are initialized at the initial weight, Wini, if set. If Weight Initial has not been set, this will be the maximum weight, Wmax.
MAXIMUM
All weights are initialized at the maximum weight, Wmax.
RANDOMIZED
Weights are uniformly distributed in Wmin,Wmax or Wini,Wmax if Weight Initial has been set.
Weight MaximumrDefault =1.0
The maximum particle weight, Wmax.
Constraint: 1.0rWmin (If Wini has been set then 1.0rWini.)
Weight MinimumrDefault =0.1
The minimum achievable weight of any particle, Wmin. Once achieved, no further weight reduction is possible.
Constraint: 0.0rWmax (If Wini has been set then 0.0rWini.)
Weight ResetaDefault =MAXIMUM
Determines how particle weights are re-initialized.
INITIAL
Weights are re-initialized at the initial weight if set. If Weight Initial has not been set, this will be the maximum weight.
MAXIMUM
Weights are re-initialized at the maximum weight.
RANDOMIZED
Weights are uniformly distributed in Wmin,Wmax or Wini,Wmax if Weight Initial has been set.
Weight ValuerDefault =0.01
The constant Wval used with Weight Decrease=INTEREST.
Constraint: 0.0r13.

12.2 Description of the SMP optional parameters

This section details additional options available to users of multithreaded implementations of the NAG Library. In particular it includes the option SMP Callback Thread Safe, which must be set before calling e05sbf with multiple threads.
SMP Callback Thread SafeaDefault =WARNING
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.
YES
The callback routines have been programmed in a thread safe way. The algorithm will use OMP_NUM_THREADS threads.
NO
The callback routines are not thread safe. Setting this option will force the algorithm to run on a single thread only, and is advisable only for debugging purposes, or if you wish to parallelize your callback functions.
WARNING
This will cause an immediate exit from e05sbf with ifail=51 if multiple threads are detected. This is to inform you that you have not declared the callback functions either to be thread safe, or that they are thread unsafe and you wish the algorithm to run in serial.
SMP Local Minimizer ExternalaDefault =ALL
Determines how many threads will attempt to locally minimize the best found solution after the routine has exited the main loop.
MASTER
Only the master thread will attempt to find any improvement. The local minimization will be launched from the best known solution. All other threads will remain effectively idle.
ALL
The master thread will perform a local minimization from the best known solution, while all other threads will perform a local minimization from randomly generated perturbations of the best known solution, increasing the chance of an improvement. Assuming all local minimizations will take approximately the same amount of computation, this will be effectively free in terms of real time. It will however increase the number of function evaluations performed.
SMP MonitoraDefault =SINGLE
SMP Monmoda
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.
SINGLE
Only one thread will invoke monmod, after all threads have performed at least one sub-iteration.
ALL
Each thread will invoke monmod each time it completes a sub-iteration. If you wish to alter x using monmod you should use this option, as monmod will only receive the arrays x, xbest, fbest and cbest private to the calling thread.
SMP SubswarmiDefault =1
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×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 i0, then this will be disabled on all the threads.
If iOMP_NUM_THREADS, then all the threads will support a particle subswarm.
SMP Thread OverruniDefault =Imax
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, the algorithm will force a synchronization between threads at the end of each iteration.