1/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2
3/*
4Copyright (C) 2015 Andres Hernandez
5
6This file is part of QuantLib, a free-software/open-source library
7for financial quantitative analysts and developers - http://quantlib.org/
8
9QuantLib is free software: you can redistribute it and/or modify it
10under the terms of the QuantLib license. You should have received a
11copy of the license along with this program; if not, please email
12<quantlib-dev@lists.sf.net>. The license is also available online at
13<http://quantlib.org/license.shtml>.
14
15This program is distributed in the hope that it will be useful, but WITHOUT
16ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
17FOR A PARTICULAR PURPOSE. See the license for more details.
18*/
19
20/*! \file particleswarmoptimization.hpp
21\brief Implementation based on:
22Clerc, M., Kennedy, J. (2002) The particle swarm-explosion, stability and
23convergence in a multidimensional complex space. IEEE Transactions on Evolutionary
24Computation, 6(2): 58–73.
25*/
26
27#ifndef quantlib_optimization_particleswarmoptimization_hpp
28#define quantlib_optimization_particleswarmoptimization_hpp
29
30#include <ql/math/optimization/problem.hpp>
31#include <ql/math/optimization/constraint.hpp>
32#include <ql/math/randomnumbers/mt19937uniformrng.hpp>
33#include <ql/experimental/math/isotropicrandomwalk.hpp>
34#include <ql/experimental/math/levyflightdistribution.hpp>
35#include <ql/math/randomnumbers/seedgenerator.hpp>
36
37#include <random>
38
39namespace QuantLib {
40
41 /*! The process is as follows:
42 M individuals are used to explore the N-dimensional parameter space:
43 \f$ X_{i}^k = (X_{i, 1}^k, X_{i, 2}^k, \ldots, X_{i, N}^k) \f$ is the kth-iteration for the ith-individual.
44
45 X is updated via the rule
46 \f[
47 X_{i, j}^{k+1} = X_{i, j}^k + V_{i, j}^{k+1}
48 \f]
49 with V being the "velocity" that updates the position:
50 \f[
51 V_{i, j}^{k+1} = \chi\left(V_{i, j}^k + c_1 r_{i, j}^k (P_{i, j}^k - X_{i, j}^k)
52 + c_2 R_{i, j}^k (G_{i, j}^k - X_{i, j}^k)\right)
53 \f]
54 where c are constants, r and R are uniformly distributed random numbers in the range [0, 1], and
55 \f$ P_{i, j} \f$ is the personal best parameter set for individual i up to iteration k
56 \f$ G_{i, j} \f$ is the global best parameter set for the swarm up to iteration k.
57 \f$ c_1 \f$ is the self recognition coefficient
58 \f$ c_2 \f$ is the social recognition coefficient
59
60 This version is known as the PSO with constriction factor (PSO-Co).
61 PSO with inertia factor (PSO-In) updates the velocity according to:
62 \f[
63 V_{i, j}^{k+1} = \omega V_{i, j}^k + \hat{c}_1 r_{i, j}^k (P_{i, j}^k - X_{i, j}^k)
64 + \hat{c}_2 R_{i, j}^k (G_{i, j}^k - X_{i, j}^k)
65 \f]
66 and is accessible from PSO-Co by setting \f$ \omega = \chi \f$,
67 and \f$ \hat{c}_{1,2} = \chi c_{1,2} \f$.
68
69 These two versions of PSO are normally referred to as canonical PSO.
70
71 Convergence of PSO-Co is improved if \f$ \chi \f$ is chosen as
72 \f$ \chi = \frac{2}{\vert 2-\phi-\sqrt{\phi^2 - 4\phi}\vert} \f$,
73 with \f$ \phi = c_1 + c_2 \f$.
74 Stable convergence is achieved if \f$ \phi >= 4 \f$. Clerc and Kennedy recommend
75 \f$ c_1 = c_2 = 2.05 \f$ and \f$ \phi = 4.1 \f$.
76
77 Different topologies can be chosen for G, e.g. instead of it being the best
78 of the swarm, it is the best of the nearest neighbours, or some other form.
79
80 In the canonical PSO, the inertia function is trivial. It is simply a
81 constant (the inertia) multiplying the previous iteration's velocity. The
82 value of the inertia constant determines the weight of a global search over
83 local search. Like in the case of the topology, other possibilities for the
84 inertia function are also possible, e.g. a function that interpolates between a
85 high inertia at the beginning of the optimization (hence prioritizing a global
86 search) and a low inertia towards the end of the optimization (hence prioritizing
87 a local search).
88
89 The optimization stops either because the number of iterations has been reached
90 or because the stationary function value limit has been reached.
91 */
92 class ParticleSwarmOptimization : public OptimizationMethod {
93 public:
94 class Inertia;
95 class Topology;
96 ParticleSwarmOptimization(Size M,
97 ext::shared_ptr<Topology> topology,
98 ext::shared_ptr<Inertia> inertia,
99 Real c1 = 2.05,
100 Real c2 = 2.05,
101 unsigned long seed = SeedGenerator::instance().get());
102 explicit ParticleSwarmOptimization(Size M,
103 ext::shared_ptr<Topology> topology,
104 ext::shared_ptr<Inertia> inertia,
105 Real omega,
106 Real c1,
107 Real c2,
108 unsigned long seed = SeedGenerator::instance().get());
109 void startState(Problem &P, const EndCriteria &endCriteria);
110 EndCriteria::Type minimize(Problem& P, const EndCriteria& endCriteria) override;
111
112 protected:
113 std::vector<Array> X_, V_, pBX_, gBX_;
114 Array pBF_, gBF_;
115 Array lX_, uX_;
116 Size M_, N_;
117 Real c0_, c1_, c2_;
118 MersenneTwisterUniformRng rng_;
119 ext::shared_ptr<Topology> topology_;
120 ext::shared_ptr<Inertia> inertia_;
121 };
122
123 //! Base inertia class used to alter the PSO state
124 /*! This pure virtual base class provides the access to the PSO state
125 which the particular inertia algorithm will change upon each iteration.
126 */
127 class ParticleSwarmOptimization::Inertia {
128 friend class ParticleSwarmOptimization;
129 public:
130 virtual ~Inertia() = default;
131 //! initialize state for current problem
132 virtual void setSize(Size M, Size N, Real c0, const EndCriteria &endCriteria) = 0;
133 //! produce changes to PSO state for current iteration
134 virtual void setValues() = 0;
135 protected:
136 ParticleSwarmOptimization *pso_;
137 std::vector<Array> *X_, *V_, *pBX_, *gBX_;
138 Array *pBF_, *gBF_;
139 Array *lX_, *uX_;
140
141 virtual void init(ParticleSwarmOptimization *pso) {
142 pso_ = pso;
143 X_ = &pso_->X_;
144 V_ = &pso_->V_;
145 pBX_ = &pso_->pBX_;
146 gBX_ = &pso_->gBX_;
147 pBF_ = &pso_->pBF_;
148 gBF_ = &pso_->gBF_;
149 lX_ = &pso_->lX_;
150 uX_ = &pso_->uX_;
151 }
152 };
153
154 //! Trivial Inertia
155 /* Inertia is a static value
156 */
157 class TrivialInertia : public ParticleSwarmOptimization::Inertia {
158 public:
159 inline void setSize(Size M, Size N, Real c0, const EndCriteria& endCriteria) override {
160 c0_ = c0;
161 M_ = M;
162 }
163 inline void setValues() override {
164 for (Size i = 0; i < M_; i++) {
165 (*V_)[i] *= c0_;
166 }
167 }
168
169 private:
170 Real c0_;
171 Size M_;
172 };
173
174 //! Simple Random Inertia
175 /* Inertia value gets multiplied with a random number
176 between (threshold, 1)
177 */
178 class SimpleRandomInertia : public ParticleSwarmOptimization::Inertia {
179 public:
180 SimpleRandomInertia(Real threshold = 0.5, unsigned long seed = SeedGenerator::instance().get())
181 : threshold_(threshold), rng_(seed) {
182 QL_REQUIRE(threshold_ >= 0.0 && threshold_ < 1.0, "Threshold must be a Real in [0, 1)");
183 }
184 inline void setSize(Size M, Size N, Real c0, const EndCriteria& endCriteria) override {
185 M_ = M;
186 c0_ = c0;
187 }
188 inline void setValues() override {
189 for (Size i = 0; i < M_; i++) {
190 Real val = c0_*(threshold_ + (1.0 - threshold_)*rng_.nextReal());
191 (*V_)[i] *= val;
192 }
193 }
194
195 private:
196 Real c0_, threshold_;
197 Size M_;
198 MersenneTwisterUniformRng rng_;
199 };
200
201 //! Decreasing Inertia
202 /* Inertia value gets decreased every iteration until it reaches
203 a value of threshold when iteration reaches the maximum level
204 */
205 class DecreasingInertia : public ParticleSwarmOptimization::Inertia {
206 public:
207 DecreasingInertia(Real threshold = 0.5)
208 : threshold_(threshold) {
209 QL_REQUIRE(threshold_ >= 0.0 && threshold_ < 1.0, "Threshold must be a Real in [0, 1)");
210 }
211 inline void setSize(Size M, Size N, Real c0, const EndCriteria& endCriteria) override {
212 N_ = N;
213 c0_ = c0;
214 iteration_ = 0;
215 maxIterations_ = endCriteria.maxIterations();
216 }
217 inline void setValues() override {
218 Real c0 = c0_*(threshold_ + (1.0 - threshold_)*(maxIterations_ - iteration_) / maxIterations_);
219 for (Size i = 0; i < M_; i++) {
220 (*V_)[i] *= c0;
221 }
222 }
223
224 private:
225 Real c0_, threshold_;
226 Size M_, N_, maxIterations_, iteration_;
227 };
228
229 //! AdaptiveInertia
230 /* Alen Lukic, Approximating Kinetic Parameters Using Particle
231 Swarm Optimization.
232 */
233 class AdaptiveInertia : public ParticleSwarmOptimization::Inertia {
234 public:
235 AdaptiveInertia(Real minInertia, Real maxInertia, Size sh = 5, Size sl = 2)
236 :minInertia_(minInertia), maxInertia_(maxInertia),
237 sh_(sh), sl_(sl) {};
238 inline void setSize(Size M, Size N, Real c0, const EndCriteria& endCriteria) override {
239 M_ = M;
240 c0_ = c0;
241 adaptiveCounter = 0;
242 best_ = QL_MAX_REAL;
243 started_ = false;
244 }
245 void setValues() override;
246
247 private:
248 Real c0_, best_;
249 Real minInertia_, maxInertia_;
250 Size M_;
251 Size sh_, sl_;
252 Size adaptiveCounter;
253 bool started_;
254 };
255
256 //! Levy Flight Inertia
257 /* As long as the particle keeps getting frequent updates to its
258 personal best value, the inertia behaves like a SimpleRandomInertia,
259 but after a number of iterations without improvement, the behaviour
260 changes to that of a Levy flight ~ u^{-1/\alpha}
261 */
262 class LevyFlightInertia : public ParticleSwarmOptimization::Inertia {
263 public:
264 LevyFlightInertia(Real alpha, Size threshold,
265 unsigned long seed = SeedGenerator::instance().get())
266 :rng_(seed), generator_(seed), flight_(generator_, LevyFlightDistribution(1.0, alpha),
267 1, Array(1, 1.0), seed),
268 threshold_(threshold) {};
269 inline void setSize(Size M, Size N, Real c0, const EndCriteria& endCriteria) override {
270 M_ = M;
271 N_ = N;
272 c0_ = c0;
273 adaptiveCounter_ = std::vector<Size>(M_, 0);
274 }
275 inline void setValues() override {
276 for (Size i = 0; i < M_; i++) {
277 if ((*pBF_)[i] < personalBestF_[i]) {
278 personalBestF_[i] = (*pBF_)[i];
279 adaptiveCounter_[i] = 0;
280 }
281 else {
282 adaptiveCounter_[i]++;
283 }
284 if (adaptiveCounter_[i] <= threshold_) {
285 //Simple Random Inertia
286 (*V_)[i] *= c0_*(0.5 + 0.5*rng_.nextReal());
287 }
288 else {
289 //If particle has not found a new personal best after threshold_ iterations
290 //then trigger a Levy flight pattern for the speed
291 flight_.nextReal<Real *>(first: &(*V_)[i][0]);
292 }
293 }
294 }
295
296 protected:
297 void init(ParticleSwarmOptimization* pso) override {
298 ParticleSwarmOptimization::Inertia::init(pso);
299 personalBestF_ = *pBF_;
300 flight_.setDimension(dim: N_, lowerBound: *lX_, upperBound: *uX_);
301 }
302
303 private:
304 MersenneTwisterUniformRng rng_;
305 std::mt19937 generator_;
306 IsotropicRandomWalk<LevyFlightDistribution, std::mt19937> flight_;
307 Array personalBestF_;
308 std::vector<Size> adaptiveCounter_;
309 Real c0_;
310 Size M_, N_;
311 Size threshold_;
312 };
313
314 //! Base topology class used to determine the personal and global best
315 /*! This pure virtual base class provides the access to the PSO state
316 which the particular topology algorithm will change upon each iteration.
317 */
318 class ParticleSwarmOptimization::Topology {
319 friend class ParticleSwarmOptimization;
320 public:
321 virtual ~Topology() = default;
322 //! initialize state for current problem
323 virtual void setSize(Size M) = 0;
324 //! produce changes to PSO state for current iteration
325 virtual void findSocialBest() = 0;
326 protected:
327 ParticleSwarmOptimization *pso_;
328 std::vector<Array> *X_, *V_, *pBX_, *gBX_;
329 Array *pBF_, *gBF_;
330 private:
331 void init(ParticleSwarmOptimization *pso) {
332 pso_ = pso;
333 X_ = &pso_->X_;
334 V_ = &pso_->V_;
335 pBX_ = &pso_->pBX_;
336 gBX_ = &pso_->gBX_;
337 pBF_ = &pso_->pBF_;
338 gBF_ = &pso_->gBF_;
339 }
340 };
341
342 //! Global Topology
343 /* The global best as seen by each particle is the best from amongst
344 all particles
345 */
346 class GlobalTopology : public ParticleSwarmOptimization::Topology {
347 public:
348 inline void setSize(Size M) override { M_ = M; }
349 inline void findSocialBest() override {
350 Real bestF = (*pBF_)[0];
351 Size bestP = 0;
352 for (Size i = 1; i < M_; i++) {
353 if (bestF < (*pBF_)[i]) {
354 bestF = (*pBF_)[i];
355 bestP = i;
356 }
357 }
358 Array& x = (*pBX_)[bestP];
359 for (Size i = 0; i < M_; i++) {
360 if (i != bestP) {
361 (*gBX_)[i] = x;
362 (*gBF_)[i] = bestF;
363 }
364 }
365 }
366
367 private:
368 Size M_;
369 };
370
371 //! K-Neighbor Topology
372 /* The global best as seen by each particle is the best from amongst
373 the previous K and next K neighbors. For particle I, the best is
374 then taken from amongst the [I - K, I + K] particles.
375 */
376 class KNeighbors : public ParticleSwarmOptimization::Topology {
377 public:
378 KNeighbors(Size K = 1) :K_(K) {
379 QL_REQUIRE(K > 0, "Neighbors need to be larger than 0");
380 }
381 inline void setSize(Size M) override {
382 M_ = M;
383 QL_ENSURE(K_ < M, "Number of neighbors need to be smaller than total particles in swarm");
384 }
385 void findSocialBest() override;
386
387 private:
388 Size K_, M_;
389 };
390
391 //! Clubs Topology
392 /* H.M. Emara, Adaptive Clubs-based Particle Swarm Optimization
393 Each particle is originally assigned to a default number of clubs
394 from among the total set. The best as seen by each particle is the
395 best from amongst the clubs to which the particle belongs.
396 Underperforming particles join more clubs randomly (up to a maximum
397 number) to widen the particles that influence them, while
398 overperforming particles leave clubs randomly (down to a minimum
399 number) to avoid early convergence to local minima.
400 */
401 class ClubsTopology : public ParticleSwarmOptimization::Topology {
402 public:
403 ClubsTopology(Size defaultClubs, Size totalClubs,
404 Size maxClubs, Size minClubs,
405 Size resetIteration, unsigned long seed = SeedGenerator::instance().get());
406 void setSize(Size M) override;
407 void findSocialBest() override;
408
409 private:
410 Size totalClubs_, maxClubs_, minClubs_, defaultClubs_;
411 Size iteration_ = 0, resetIteration_;
412 Size M_;
413 std::vector<std::vector<bool> > clubs4particles_;
414 std::vector<std::vector<bool> > particles4clubs_;
415 std::vector<Size> bestByClub_;
416 std::vector<Size> worstByClub_;
417 std::mt19937 generator_;
418 std::uniform_int_distribution<QuantLib::Size> distribution_;
419 using param_type = decltype(distribution_)::param_type;
420
421 void leaveRandomClub(Size particle, Size currentClubs);
422 void joinRandomClub(Size particle, Size currentClubs);
423 };
424
425}
426
427#endif
428

source code of quantlib/ql/experimental/math/particleswarmoptimization.hpp