Particle Swarm Optimization (PSO) is a stochastic swarm intelligence algorithm for global optimization of potentially ill-behaved nonlinear functions.


Definition and Syntax

bool pso(arma::vec& init_out_vals, std::function<double (const arma::vec& vals_inp, arma::vec* grad_out, void* opt_data)> opt_objfn, void* opt_data);
bool pso(arma::vec& init_out_vals, std::function<double (const arma::vec& vals_inp, arma::vec* grad_out, void* opt_data)> opt_objfn, void* opt_data, algo_settings_t& settings);

Function arguments:

  • init_out_vals a column vector of initial values; will be replaced by the solution values.
  • opt_objfn the function to be minimized, taking three arguments:
    • vals_inp a vector of inputs;
    • grad_out an empty vector, as PSO does not require the gradient to be known/exist; and
    • opt_data additional parameters passed to the function.
  • opt_data additional parameters passed to the function.
  • settings parameters controlling the optimization routine; see below.

Optimization control parameters:

  • bool vals_bound whether the search space is bounded. If true, then
    • arma::vec lower_bounds this defines the lower bounds.
    • arma::vec upper_bounds this defines the upper bounds.

  • int pso_n_pop population size of each generation.
  • int pso_n_gen number of vectors to generate.
  • int pso_check_freq number of generations between convergence checks.

  • int pso_inertia_method method of inertia decay:
    • pso_inertia_method = 1 linear decay between pso_par_w_max and pso_par_w_min.
    • pso_inertia_method = 2 dampening using pso_par_w_damp parameter.

  • int pso_velocity_method method of updating the velocity weights:
    • pso_velocity_method = 1 fixed weights pso_par_c_cog and pso_par_c_soc.
    • pso_velocity_method = 2 linear decay between pso_par_initial_c_cog and pso_par_final_c_cog, and initial_c_soc and pso_par_final_c_soc.
  • double pso_par_c_cog weight value on the 'cognitive' factor.
  • double pso_par_c_soc weight value on the 'social' factor.

  • arma::vec pso_initial_lb lower bounds on the initial population; defaults to init_out_vals $- \ 0.5$.
  • arma::vec pso_initial_ub upper bounds on the initial population; defaults to init_out_vals $+ \ 0.5$.

Details

The PSO method in OptimLib is available in two varieties: the simple version, outlined first, and the a hybrid DE-PSO algorithm.


Basic Particle Swarm Optimization (PSO)

Let $X^{(i)}$ denote a $N_p \times d$-dimensional array of values at stage $i$ of the algorithm. The basic PSO algorithm is as follows.

  • Update the velocity and position Sample $R_C,R_S \sim U^d(0,1)$, $d$-dimensional iid uniform random vectors, and set $$V^{(i+1)}(j.:) = w V^{(i+1)}(j,:) + c_C R_C \odot (X_b^{(i)} (j,:) - X^{(i)}(j,:)) + c_S R_S \odot (g_b - X^{(i)}(j,:))$$ $$X^{(i+1)}(j,:) = X^{(i)}(j,:) + V^{(i+1)}(j,:)$$ where the values $c_C$ and $c_S$ are set by pso_par_c_cog and pso_par_c_soc, respectively, and '$\odot$' denotes the element-by-element (Hadamard) product.
  • Update local-best particle
  • $$X_b^{(i+1)} (j,:) = \begin{cases} X^{(i+1)}(j,:) & \text{ if } f(X^{(i+1)}(j,:)) < f(X_b^{(i)}) \\ X_b^{(i)} (j,:) & \text{ else } \end{cases}$$
  • Update the global-best particle
  • Let $$X^* := \arg \min_{X_j} f(X(j,:))$$ Then $$g_b = \begin{cases} X^* & \text{ if } f(X^*) < f(g_b) \\ g_b & \text{ else } \end{cases}$$
  • Reduce $w$ and update weights

The algorithm stops when the relative improvement in the objective function is less than err_tol between pso_check_freq number of generations, or when the total number of 'generations' exceeds n_gen.


Particle Swarm Optimization with Differentially-Perturbed Velocity (PSO-DV)



Examples

To illustrate the effectiveness of PSO, we will tackle two well-known performance tests from the numerical optimization literature:

  • The Table function:
  • $$\min_{x \in [-10,10]^2} \left\{ - \left| \sin (x_1) \cos(x_2) \exp \left( \left| 1 - \frac{\sqrt{x_1^2 + x_2^2}}{\pi} \right| \right) \right| \right\}$$
  • Bukin's function No.6:
  • $$\min_{x \in [-15,-5], y \in [-3,3]} \left\{ 100 \sqrt{| y - 0.01 x^2|} + 0.01 | x + 10 | \right\}$$

The first function is bumpy, contains many local minima, and four global minima; the second contains a particularly difficult-to-traverse ridge. Plots of both functions are given below.

Table Function

Bukin No. 6


Code implementing these examples with OptimLib is given below.

#include "optim.hpp"

double table_fn(const arma::vec& vals_inp, arma::vec* grad_out, void* opt_data)
{
    const double x = vals_inp(0);
    const double y = vals_inp(1);
    const double pi = arma::datum::pi;

    double obj_val = - std::abs( std::sin(x)*std::cos(y)*std::exp( std::abs(1.0 - std::sqrt(x*x + y*y)/pi) ) );
    //
    return obj_val;
}

double bukin_fn(const arma::vec& vals_inp, arma::vec* grad_out, void* opt_data)
{
    const double x = vals_inp(0);
    const double y = vals_inp(1);

    double obj_val = 100*std::sqrt(std::abs(y - 0.01*x*x)) + 0.01*std::abs(x + 10);
    //
    return obj_val;
}

int main()
{
    //
    // Table function

    optim::algo_settings_t settings_1;

    settings_1.pso_center_particle = false;
    settings_1.pso_par_bounds = true;

    arma::vec x_1 = arma::zeros(2,1);

    settings_1.pso_lb = arma::zeros(2,1) - 10.0;
    settings_1.pso_ub = arma::zeros(2,1) + 10.0;

    settings_1.pso_n_pop = 5000;
    settings_1.pso_n_gen = 4000;

    bool success_1 = optim::pso(x_1,table_fn,nullptr,settings_1);

    if (success_1) {
        std::cout << "pso: Table test completed successfully." << std::endl;
    } else {
        std::cout << "pso: Table test completed unsuccessfully." << std::endl;
    }

    arma::cout << "pso: solution to Table test:\n" << x_1 << arma::endl;

    //
    // Bukin function no. 6

    arma::vec x_2 = arma::ones(2,1) + 1.0; // initial values: (2,2)

    bool success_2 = optim::de(x_2,bukin_fn,nullptr);

    if (success_2) {
        std::cout << "de: Bukin test completed successfully." << std::endl;
    } else {
        std::cout << "de: Bukin test completed unsuccessfully." << std::endl;
    }

    arma::cout << "de: solution to Bukin test:\n" << x_2 << arma::endl;

    return 0;
}