31 #include <boost/math/special_functions/binomial.hpp>
32 #include <boost/random/uniform_real.hpp>
33 #include <boost/random/uniform_int.hpp>
34 #include <boost/random/variate_generator.hpp>
36 #include "../exceptions.h"
37 #include "../population.h"
38 #include "../problem/base.h"
39 #include "../population.h"
40 #include "../util/neighbourhood.h"
45 namespace pagmo {
namespace algorithm {
62 nspso::nspso(
int gen,
double minW,
double maxW,
double C1,
double C2,
63 double CHI,
double v_coeff,
int leader_selection_range,
72 m_leader_selection_range(leader_selection_range),
73 m_diversity_mechanism(diversity_mechanism)
76 pagmo_throw(value_error,
"number of generations must be nonnegative");
79 if (minW <= 0 || maxW <= 0 || C1 <=0 || C2 <= 0 || CHI <=0) {
80 pagmo_throw(value_error,
"minW, maxW, C1, C2 and CHI should be greater than 0");
84 pagmo_throw(value_error,
"minW should be lower than maxW");
87 if (v_coeff <= 0 || v_coeff > 1) {
88 pagmo_throw(value_error,
"v_coeff should be in the ]0,1] range");
91 if (leader_selection_range <=0 || leader_selection_range > 100) {
92 pagmo_throw(value_error,
"leader_selection_range should be in the ]0,100] range");
120 pagmo_throw(value_error,
"There is no continuous part in the problem decision vector for NSPSO to optimise");
123 if( prob_c_dimension != 0 ){
124 pagmo_throw(value_error,
"The problem is not box constrained and NSPSO is not suitable to solve it");
127 if( prob_f_dimension < 2 ){
128 pagmo_throw(value_error,
"The problem is not multi-objective. Use a single-objectie optimization algorithm instead");
132 pagmo_throw(value_error,
"non existing diversity mechanism method");
136 if (NP == 0 || m_gen == 0) {
144 double v_width = (ub[d] - lb[d]) * m_v_coeff;
145 minV[d] = -1.0 * v_width;
150 std::vector<nspso_individual> nextPopList;
151 for(
unsigned int i=0; i < NP; ++i) {
152 nextPopList.push_back(nspso_individual());
162 for(
int g = 0; g < m_gen; ++g) {
164 std::vector<population::size_type> bestNonDomIndices;
165 std::vector<fitness_vector> fit(NP);
166 std::vector<constraint_vector> cons(NP);
168 fit[i] = nextPopList[i].best_f;
169 cons[i] = nextPopList[i].best_c;
174 std::vector<std::vector<population::size_type> > dom_list = compute_domination_list(prob,fit,cons);
175 std::vector<population::size_type> pareto_rank = compute_pareto_rank(dom_list);
176 std::vector<std::vector<population::size_type> > pareto_fronts = compute_pareto_fronts(pareto_rank);
177 std::vector<double> crowding_d = compute_crowding_d(fit, pareto_fronts);
179 crowding_pareto_comp comp(pareto_rank, crowding_d);
180 std::vector<population::size_type> dummy(NP);
181 for(
unsigned int i=0; i<NP; ++i) dummy[i] = i;
182 std::sort(dummy.begin(), dummy.end(),comp);
183 if(pareto_fronts[0].size() > 1) {
184 bestNonDomIndices = std::vector<population::size_type>(dummy.begin(), dummy.begin() + pareto_fronts[0].size());
186 bestNonDomIndices = std::vector<population::size_type>(dummy.begin(), dummy.begin() + 2);
190 std::vector<decision_vector> nonDomChromosomes;
192 std::vector<std::vector<population::size_type> > dom_list = compute_domination_list(prob,fit,cons);
193 std::vector<population::size_type> pareto_rank = compute_pareto_rank(dom_list);
194 std::vector<std::vector<population::size_type> > pareto_fronts = compute_pareto_fronts(pareto_rank);
196 for(
unsigned int i = 0; i < pareto_fronts[0].size(); ++i) {
197 nonDomChromosomes.push_back(nextPopList[pareto_fronts[0][i]].best_x);
200 std::vector<double> nadir = compute_nadir(fit, pareto_rank);
201 std::vector<double> ideal = compute_ideal(fit, pareto_rank);
205 if (prob_f_dimension == 2) {
206 delta = ( (nadir[0] - ideal[0]) + (nadir[1] - ideal[1]) ) / (nonDomChromosomes.size()-1);
207 }
else if (prob_f_dimension == 3) {
208 const double d1 = nadir[0] - ideal[0];
209 const double d2 = nadir[1] - ideal[1];
210 const double d3 = nadir[2] - ideal[2];
211 const double N = nonDomChromosomes.size();
212 delta = sqrt(4*d2*d1*N + 4*d3*d1*N + 4*d2*d3*N + pow(d1,2) + pow(d2,2) + pow(d3,2)
213 - 2*d2*d1 - 2*d3*d1 - 2*d2*d3 + d1 + d2 + d3) / (2*(N-1));
215 for(
unsigned int i=0; i<nadir.size(); ++i) {
216 delta *= nadir[i] - ideal[i];
218 delta = pow(delta, 1.0/nadir.size())/nonDomChromosomes.size();
221 std::vector<int> count(nonDomChromosomes.size(),0);
222 compute_niche_count(count, nonDomChromosomes, delta);
225 if(pareto_fronts[0].size() > 1) {
226 for(
unsigned int i=0; i < tmp.size(); ++i) {
227 bestNonDomIndices.push_back(pareto_fronts[0][tmp[i]]);
230 unsigned int minPopSize = 2;
231 for(
unsigned f = 0; minPopSize > 0 && f < pareto_fronts.size(); ++f) {
232 for(
unsigned int i=0; minPopSize > 0 && i < pareto_fronts[f].size(); ++i) {
233 bestNonDomIndices.push_back(pareto_fronts[f][i]);
239 std::vector<double> maxmin(NP,0);
240 compute_maxmin(maxmin, fit);
243 for(;i < bestNonDomIndices.size() && maxmin[bestNonDomIndices[i]] < 0; ++i);
245 bestNonDomIndices = std::vector<population::size_type>(bestNonDomIndices.begin(), bestNonDomIndices.begin() + i);
249 const double W = m_maxW - (m_maxW-m_minW)/m_gen * g;
258 int ext = ceil(bestNonDomIndices.size()*m_leader_selection_range/100.0)-1;
259 if (ext < 1) ext = 1;
260 unsigned int leaderIdx;
262 leaderIdx = boost::uniform_int<int>(0,ext)(
m_drng);
263 }
while (bestNonDomIndices[leaderIdx] == idx);
264 decision_vector leader = nextPopList[bestNonDomIndices[leaderIdx]].best_x;
267 const double r1 = boost::uniform_real<double>(0,1)(
m_drng);
268 const double r2 = boost::uniform_real<double>(0,1)(
m_drng);
273 for(decision_vector::size_type i = 0; i < cur_X.size(); ++i) {
274 double v = W*cur_V[i] +
275 m_C1*r1*(best_X[i] - cur_X[i]) +
276 m_C2*r2*(leader[i] - cur_X[i]);
281 else if(v < minV[i]) {
285 double x = cur_X[i] + m_CHI*v;
289 }
else if (x < lb[i]) {
299 nextPopList.push_back(nspso_individual());
300 nextPopList[idx+NP].cur_x = newX;
301 nextPopList[idx+NP].best_x = newX;
302 nextPopList[idx+NP].cur_v = newV;
303 nextPopList[idx+NP].cur_f = prob.
objfun(newX);
304 nextPopList[idx+NP].best_f = nextPopList[idx+NP].cur_f;
306 nextPopList[idx+NP].best_c = nextPopList[idx+NP].cur_c;
310 std::vector<fitness_vector> nextPop_fit(nextPopList.size());
311 std::vector<constraint_vector> nextPop_cons(nextPopList.size());
312 for(
unsigned int i=0; i<nextPopList.size(); ++i) {
313 nextPop_fit[i] = nextPopList[i].best_f;
314 nextPop_cons[i] = nextPopList[i].best_c;
317 std::vector<population::size_type> bestNextPopIndices(NP,0);
319 if(m_diversity_mechanism !=
MAXMIN) {
320 std::vector<std::vector<population::size_type> > nextPop_pareto_fronts = compute_pareto_fronts(prob, nextPop_fit, nextPop_cons);
321 for(
unsigned int f = 0, i=0; i<NP && f < nextPop_pareto_fronts.size(); ++f) {
322 if(nextPop_pareto_fronts[f].size() < NP-i) {
323 for(
unsigned int j = 0; j < nextPop_pareto_fronts[f].size(); ++j) {
324 bestNextPopIndices[i] = nextPop_pareto_fronts[f][j];
328 boost::uniform_int<int> pop_idx(0,nextPop_pareto_fronts[f].size());
329 boost::variate_generator<boost::mt19937 &, boost::uniform_int<int> > p_idx(
m_urng,pop_idx);
330 std::random_shuffle(nextPop_pareto_fronts[f].begin(), nextPop_pareto_fronts[f].end(), p_idx);
331 for(
unsigned int j = 0; i<NP; ++j) {
332 bestNextPopIndices[i] = nextPop_pareto_fronts[f][j];
338 std::vector<double> maxmin(2*NP,0);
339 compute_maxmin(maxmin, nextPop_fit);
341 bestNextPopIndices = std::vector<population::size_type>(tmp.begin(), tmp.begin() + NP);
345 std::vector<nspso_individual> tmpPop(NP);
346 for(
unsigned int i=0; i < NP; ++i) {
347 tmpPop[i] = nextPopList[bestNextPopIndices[i]];
350 nextPopList = tmpPop;
356 pop.push_back(nextPopList[i].best_x);
357 pop.set_x(i, nextPopList[i].cur_x);
358 pop.set_v(i, nextPopList[i].cur_v);
363 double nspso::minfit(
unsigned int i,
unsigned int j,
const std::vector<fitness_vector> &fit)
const
365 double min = fit[i][0] - fit[j][0];
366 for(
unsigned int f=0; f<fit[i].size(); ++f) {
367 double tmp = fit[i][f] - fit[j][f];
375 void nspso::compute_maxmin(std::vector<double> &maxmin,
const std::vector<fitness_vector> &fit)
const
377 const unsigned int NP = fit.size();
378 for(
unsigned int i=0; i<NP; ++i) {
379 maxmin[i] = minfit(i, (i+1)%NP, fit);
380 for(
unsigned j=0; j<NP; ++j) {
382 double tmp = minfit(i, j, fit);
383 if(tmp > maxmin[i]) {
393 double nspso::euclidian_distance(
const std::vector<double> &x,
const std::vector<double> &y)
const
395 pagmo_assert(x.size() == y.size());
397 for(
unsigned int i = 0; i < x.size(); ++i) {
398 sum+= pow(x[i]-y[i],2);
403 void nspso::compute_niche_count(std::vector<int> &count,
const std::vector<std::vector<double> > &chromosomes,
double delta)
const
405 std::fill(count.begin(), count.end(),0);
406 for(
unsigned int i=0; i<chromosomes.size(); ++i) {
407 for(
unsigned int j=0; j<chromosomes.size(); ++j) {
408 if(euclidian_distance(chromosomes[i], chromosomes[j]) < delta) {
416 std::vector<std::vector<population::size_type> > nspso::compute_domination_list(
const pagmo::problem::base &prob,
417 const std::vector<fitness_vector> &fit,
418 const std::vector<constraint_vector> &cons)
const
420 std::vector<population::size_type> dummy;
421 std::vector<std::vector<population::size_type> > domination_list(fit.size(), dummy);
423 for(
unsigned int i=0; i<fit.size();++i) {
424 for(
unsigned int j=0; j<fit.size(); ++j) {
426 if(prob.
compare_fc(fit[i],cons[i],fit[j],cons[j])) {
427 domination_list[i].push_back(j);
432 return domination_list;
435 std::vector<population::size_type> nspso::compute_domination_count(
const std::vector<std::vector<population::size_type> > &dom_list)
const
437 std::vector<population::size_type> domination_count(dom_list.size(),0);
438 for(
unsigned int i=0; i<dom_list.size(); ++i) {
439 for(
unsigned int j = 0; j < dom_list[i].size(); ++j) {
440 domination_count[dom_list[i][j]]++;
444 return domination_count;
447 std::vector<population::size_type> nspso::compute_pareto_rank(
const std::vector<std::vector<population::size_type> > &dom_list)
const
449 std::vector<population::size_type> pareto_rank(dom_list.size(),0);
452 std::vector<population::size_type> F,S;
455 std::vector<population::size_type> dom_count = compute_domination_count(dom_list);
459 if (dom_count[idx] == 0) {
464 unsigned int irank = 1;
467 while (F.size()!=0) {
472 dom_count[dom_list[F[i]][j]]--;
473 if (dom_count[dom_list[F[i]][j]] == 0){
474 S.push_back(dom_list[F[i]][j]);
475 pareto_rank[dom_list[F[i]][j]] = irank;
488 std::vector<double> nspso::compute_crowding_d(
const std::vector<fitness_vector> &fit,
const std::vector<std::vector<population::size_type> > &pareto_fronts)
const {
490 std::vector<double> crowding_d(fit.size(),0);
492 for(
unsigned f=0; f < pareto_fronts.size(); ++f) {
493 std::vector<fitness_vector::size_type> I(pareto_fronts[f]);
495 fitness_vector::size_type lastidx = I.size() - 1;
498 one_dim_fit_comp funct(fit,0);
501 for (fitness_vector::size_type i = 0; i < fit[0].size(); ++i) {
504 std::sort(I.begin(),I.end(), funct);
506 crowding_d[I[0]] = std::numeric_limits<double>::max();
507 crowding_d[I[lastidx]] = std::numeric_limits<double>::max();
509 double df = fit[I[lastidx]][i] - fit[I[0]][i];
512 crowding_d[I[j]] += 0.0;
514 crowding_d[I[j]] += (fit[I[j+1]][i] -fit[I[j-1]][i])/df;
523 std::vector<std::vector<population::size_type> > nspso::compute_pareto_fronts(
const pagmo::problem::base &prob,
524 const std::vector<fitness_vector> &fit,
525 const std::vector<constraint_vector> &cons)
const
527 std::vector<std::vector<population::size_type> > dom_list = compute_domination_list(prob,fit,cons);
528 std::vector<population::size_type> pareto_rank = compute_pareto_rank(dom_list);
530 return compute_pareto_fronts(pareto_rank);
533 std::vector<std::vector<population::size_type> > nspso::compute_pareto_fronts(
const std::vector<population::size_type> &pareto_rank)
const
535 std::vector<std::vector<population::size_type> > retval;
538 if (pareto_rank[idx] >= retval.size()) {
539 retval.resize(pareto_rank[idx] + 1);
541 retval[pareto_rank[idx]].push_back(idx);
547 fitness_vector nspso::compute_ideal(
const std::vector<fitness_vector> &fit,
const std::vector<population::size_type> &pareto_rank)
const {
549 unsigned int firstFrontIdx = 0;
550 for(;pareto_rank[firstFrontIdx] > 0; ++firstFrontIdx);
554 if (pareto_rank[idx] == 0) {
555 for(fitness_vector::size_type i = 0; i < ideal.size(); ++i) {
556 if (fit[idx][i] < ideal[i]) {
557 ideal[i] = fit[idx][i];
565 fitness_vector nspso::compute_nadir(
const std::vector<fitness_vector> &fit,
const std::vector<population::size_type> &pareto_rank)
const {
566 unsigned int firstFrontIdx = 0;
567 for(;pareto_rank[firstFrontIdx] > 0; ++firstFrontIdx);
571 if (pareto_rank[idx] == 0) {
572 for(fitness_vector::size_type i = 0; i < nadir.size(); ++i) {
573 if (fit[idx][i] > nadir[i]) {
574 nadir[i] = fit[idx][i];
586 return "Non-dominated Sorting Particle Swarm Optimizer (NSPSO)";
595 std::ostringstream s;
596 s <<
"gen:" << m_gen <<
' ';
597 s <<
"minW:" << m_minW <<
' ';
598 s <<
"maxW:" << m_maxW <<
' ';
599 s <<
"C1:" << m_C1 <<
' ';
600 s <<
"C2:" << m_C2 <<
' ';
601 s <<
"CHI:" << m_CHI <<
' ';
602 s <<
"v_coeff:" << m_v_coeff <<
' ';
603 s <<
"leader_selection_range:" << m_leader_selection_range <<
' ';
604 s <<
"diversity_mechanism: ";
605 switch (m_diversity_mechanism)
611 case MAXMIN : s <<
"MAXMIN" <<
' ';
boost::shared_ptr< base > base_ptr
Alias for shared pointer to base algorithm.
std::vector< double > decision_vector
Decision vector type.
fitness_vector cur_f
Current fitness vector.
const individual_type & get_individual(const size_type &) const
Get constant reference to individual at position n.
constraint_vector cur_c
Current constraint vector.
base_ptr clone() const
Clone method.
std::vector< population::size_type > order(const std::vector< T > &values)
Sort according the the values in the values vector but return the permutation.
size_type get_dimension() const
Return global dimension.
fitness_vector objfun(const decision_vector &) const
Return fitness of pagmo::decision_vector.
std::string human_readable_extra() const
Extra human readable algorithm info.
size_type get_i_dimension() const
Return integer dimension.
void evolve(population &) const
Evolve implementation.
nspso(int gen=100, double minW=0.4, double maxW=1.0, double C1=2.0, double C2=2.0, double CHI=1.0, double v_coeff=0.5, int leader_selection_range=10, diversity_mechanism_type=CROWDING_DISTANCE)
Constructor.
fitness_vector best_f
Best fitness vector so far.
std::vector< double > fitness_vector
Fitness vector type.
constraint_vector compute_constraints(const decision_vector &) const
Compute constraints and return constraint vector.
diversity_mechanism_type
Mechanism used to asses diversity.
c_size_type get_c_dimension() const
Return global constraints dimension.
decision_vector cur_v
Current velocity vector.
const decision_vector & get_ub() const
Upper bounds getter.
container_type::size_type size_type
Population size type.
decision_vector cur_x
Current decision vector.
rng_uint32 m_urng
Random number generator for unsigned integer values.
f_size_type get_f_dimension() const
Return fitness dimension.
constraint_vector best_c
Best constraint vector so far.
bool compare_fc(const fitness_vector &, const constraint_vector &, const fitness_vector &, const constraint_vector &) const
Simultaneous fitness-constraint comparison.
const decision_vector & get_lb() const
Lower bounds getter.
rng_double m_drng
Random number generator for double-precision floating point values.
decision_vector best_x
Best decision vector so far.
decision_vector::size_type size_type
Problem's size type: the same as pagmo::decision_vector's size type.
std::string get_name() const
Algorithm name.
Non-dominated Sorting Particle Swarm Optimizer (NSPSO)
The crowding distance is used.