# particle_swarm_optimization.py
#
# Created: Sep. 2019, M. Clarke
# ----------------------------------------------------------------------
# Imports
# ----------------------------------------------------------------------
# RCAIDE imports
import numpy as np
import scipy as sp
[docs]
def particle_swarm_optimization(func, lb, ub, ieqcons=[], f_ieqcons=None, args=(), kwargs={},
swarmsize=100, omega=0.5, phip=0.5, phig=0.5, maxiter=100,
minstep=1e-8, minfunc=1e-8, debug=False):
"""
This function perform a particle swarm optimization (PSO)
Source:
Pyswarm: https://github.com/tisimst/pyswarm
Inputs:
func : The function to be minimized [function]
lb : The lower bounds of the design variable(s) [array]
ub : The upper bounds of the design variable(s) [array]
ieqcons : A list of functions of length n such that ieqcons[j](x,*args)
>= 0.0 in a successfully optimized problem (Default: []) [list]
f_ieqcons : Returns a 1-D array in which each element must be greater or equal
to 0.0 in a successfully optimized problem. If f_ieqcons is specified,
ieqcons is ignored (Default: None) [function]
args : Additional arguments passed to objective and constraint functions [tuple]
kwargs : Additional keyword arguments passed to objective and constraint functions [dict]
swarmsize : The number of particles in the swarm (Default: 100) [int]
omega : Particle velocity scaling factor (Default: 0.5) [float]
phip : Scaling factor to search away from the particle's best known position (Default: 0.5) [scalar]
phig : Scaling factor to search away from the swarm's best known position (Default: 0.5) [scalar]
maxiter : The maximum number of iterations for the swarm to search (Default: 100) [int]
minstep : The minimum stepsize of swarm's best position before the search terminates (Default: 1e-8) [scalar]
minfunc : The minimum change of swarm's best objective value before the search terminates (Default: 1e-8) [scalar]
debug : If True, progress statements will be displayed every iteration (Default: False) [boolean]
Outputs:
g : The swarm's best known position (optimal design) [list]
f : The objective value at ``g`` [float]
Properties Used:
None
"""
# License
'''
Copyright (c) 2015, Sebastian M. Castillo Hair
All rights reserved.
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
1. Redistributions of source code must retain the above copyright
notice, this list of conditions and the following disclaimer.
2. Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in the
documentation and/or other materials provided with the distribution.
3. All advertising materials mentioning features or use of this software
must display the following acknowledgement:
This product includes software developed by Sebastian M. Castillo Hair.
4. Neither the name of Sebastian M. Castillo Hair nor the
names of its contributors may be used to endorse or promote products
derived from this software without specific prior written permission.
THIS SOFTWARE IS PROVIDED BY SEBASTIAN M. CASTILLO HAIR ''AS IS'' AND ANY
EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL SEBASTIAN M. CASTILLO HAIR BE LIABLE FOR ANY
DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.'''
assert len(lb)==len(ub), 'Lower- and upper-bounds must be the same length'
assert hasattr(func, '__call__'), 'Invalid function handle'
lb = np.array(lb)
ub = np.array(ub)
assert np.all(ub>lb), 'All upper-bound values must be greater than lower-bound values'
vhigh = np.abs(ub - lb)
vlow = -vhigh
# Check for constraint function(s) #########################################
obj = lambda x: func(x, *args, **kwargs)
if f_ieqcons is None:
if not len(ieqcons):
if debug:
print('No constraints given.')
cons = lambda x: np.array([0])
else:
if debug:
print('Converting ieqcons to a single constraint function')
cons = lambda x: np.array([y(x, *args, **kwargs) for y in ieqcons])
else:
if debug:
print('Single constraint function given in f_ieqcons')
cons = lambda x: np.array(f_ieqcons(x, *args, **kwargs))
def is_feasible(x):
check = np.all(cons(x)>=0)
return check
# Initialize the particle swarm ############################################
S = swarmsize
D = len(lb) # the number of dimensions each particle has
x = np.random.rand(S, D) # particle positions
v = np.zeros_like(x) # particle velocities
p = np.zeros_like(x) # best particle positions
fp = np.zeros(S) # best particle function values
g = [] # best swarm position
fg = 1e100 # artificial best swarm position starting value
for i in range(S):
# Initialize the particle's position
x[i, :] = lb + x[i, :]*(ub - lb)
# Initialize the particle's best known position
p[i, :] = x[i, :]
# Calculate the objective's value at the current particle's
fp[i] = obj(p[i, :])
# At the start, there may not be any feasible starting point, so just
# give it a temporary "best" point since it's likely to change
if i==0:
g = p[0, :].copy()
# If the current particle's position is better than the swarm's,
# update the best swarm position
if fp[i]<fg and is_feasible(p[i, :]):
fg = fp[i]
g = p[i, :].copy()
# Initialize the particle's velocity
v[i, :] = vlow + np.random.rand(D)*(vhigh - vlow)
# Iterate until termination criterion met ##################################
it = 1
while it<=maxiter:
rp = np.random.uniform(size=(S, D))
rg = np.random.uniform(size=(S, D))
for i in range(S):
# Update the particle's velocity
v[i, :] = omega*v[i, :] + phip*rp[i, :]*(p[i, :] - x[i, :]) + \
phig*rg[i, :]*(g - x[i, :])
# Update the particle's position, correcting lower and upper bound
# violations, then update the objective function value
x[i, :] = x[i, :] + v[i, :]
mark1 = x[i, :]<lb
mark2 = x[i, :]>ub
x[i, mark1] = lb[mark1]
x[i, mark2] = ub[mark2]
fx = obj(x[i, :])
# Compare particle's best position (if constraints are satisfied)
if fx<fp[i] and is_feasible(x[i, :]):
p[i, :] = x[i, :].copy()
fp[i] = fx
# Compare swarm's best position to current particle's position
# (Can only get here if constraints are satisfied)
if fx<fg:
if debug:
print('New best for swarm at iteration {:}: {:} {:}'.format(it, x[i, :], fx))
tmp = x[i, :].copy()
stepsize = np.sqrt(np.sum((g-tmp)**2))
if np.abs(fg - fx)<=minfunc:
print('Stopping search: Swarm best objective change less than {:}'.format(minfunc))
return tmp, fx
elif stepsize<=minstep:
print('Stopping search: Swarm best position change less than {:}'.format(minstep))
return tmp, fx
else:
g = tmp.copy()
fg = fx
if debug:
print('Best after iteration {:}: {:} {:}'.format(it, g, fg))
it += 1
print('Stopping search: maximum iterations reached --> {:}'.format(maxiter))
if not is_feasible(g):
print("However, the optimization couldn't find a feasible design. Sorry")
return g, fg
