单目标 JADE

今天来发表一篇博客吧，写的是JADE（改进的养分演化算法），是基于这篇论文实现的“http://150.214.190.154/docencia/sf1/2009-Zhang-IEEETEC.pdf”。

大概讲讲吧，写的久了，也忘的差不多了，也当作个回忆复习吧。

我全部是用STL实现的，所以不懂STL的可能就比较难看的懂了。

试函数我用的是第一个，即  

测试数据也是第一个，即1500代，30维向量，种群大小为100

测试结果和预期结果一致，精度为e-60，非常接近零了~

下面就直接上代码了

要是有不懂的或者哪里有错误的，可以直接回复或者发邮件到louzhang.swk@gmail.com

忘了说明一点了，正态分布和柯西分布那两个函数是直接用Jjj学长教的方法，用的是tr1命名空间里封装好的函数

#include<iostream>
#include<cstdio>
#include<cmath>
#include<ctime>
#include<cstdlib>
#include<algorithm>
#include <random>
#include <set>
#include <vector>
using namespace std;
 
 
 
#define PI acos(-1.0)
#define D 30                        //基因数
#define NP 100                      //种群个数
#define G 1500                      //代数
#define c 0.1
#define p 0.1
 
//http://150.214.190.154/docencia/sf1/2009-Zhang-IEEETEC.pdf
//开发环境：VS 2010 + SP1补丁
 
 
struct individual
{
    double xreal[D+10];
    double fitness;
   
    friend bool operator < (const individual& a,const individual& b)
    {
       return a.fitness-b.fitness<0;
       //individual x=a,y=b;
       //return calc_indi_fitness(x)<calc_indi_fitness(y);
    }
    //friend double calc_indi_fitness(individual& ind);
    bool individual::operator ==(const individual& other)
    {
       bool isEqul = true;
       for(int i = 0; i<D; i++)
       {
           if(this->xreal[i] != other.xreal[i])
           {
              isEqul = false;
              break;
           }
       }
       return isEqul;
    }
 
}indi;
 
 
vector<individual>P,A;
double uCR=0.5,uF=0.5;
double CR,F;
int x_best,r1,r2,xi;
 
double calc_indi_fitness(individual& ind)
{
    double sum=0.0;
    for(int i=0;i<D;i++)sum+=ind.xreal[i]*ind.xreal[i];
    return sum;
}
void init_indi(tr1::mt19937& gen)
{
    tr1::uniform_real<double>  dist(-100.0,100.0);
    for(int i=0;i<NP;i++)
    {
       for(int j = 0; j<D; j++)
       {
           indi.xreal[j] = dist(gen);
       }
       indi.fitness=calc_indi_fitness(indi);
       P.push_back(indi);
    }
}
double rand_normal(double u,double r,tr1::mt19937& gen)
{
    tr1::normal_distribution<double>  dist(u,r);
    double result = dist(gen);
    if(result > 1.0)result = 1.0;
    if(result < 0.0)result = 0.0;
    return result;
}
 
double rand_cauthy(double x0,double r,tr1::mt19937& gen)
{
    tr1::uniform_real<double>  dist(-0.5*PI,0.5*PI);
    double result = 0.0;
    while(1)
    {
       result = x0 + r * tan(dist(gen));
       if(result <= 0.0)
           continue;
       if(result >= 1.0)
       {
           result = 1.0;
           break;
       }else break;
    }
    return result;
}
int rand_int(int lb,int ub,tr1::mt19937& gen)
{
    tr1::uniform_int<int> dist(lb,ub);
    return dist(gen);
}
double rand_real(double lb,double ub,tr1::mt19937& gen)
{
    tr1::uniform_real<double>  dist(lb,ub);
    return dist(gen);
}
double Mean_A(const vector<double>& scr)
{
    double result = 0;
    vector<double>::const_iterator it = scr.begin();
    int count = 0;
    for(; it!=scr.end();++it)
    {
       result += *it;
       count++;
    }
    return (result/count);
}
 
double Mean_L(const vector<double>& sf)
{
    double sum = 0.0;
    double result = 0.0;
    double temp;
    vector<double>::const_iterator it = sf.begin();
    for(; it!=sf.end();++it)
    {
       temp = *it;
       sum += temp;
       result += temp*temp;
    }
    return (result/sum);
}
bool cmp(const individual& a,const individual& b)
{
    individual x=a,y=b;
    return calc_indi_fitness(x)<calc_indi_fitness(y);
}
int main()
{
    //int start_time=clock();
    tr1::mt19937 gen;
    gen.seed((unsigned long)time(0));
 
 
    init_indi(gen);
 
    int count=0;
    for(int g=0;g<G;g++)
    {
       //cout<<g<<endl;
       vector<double>SF,SCR;
       for (int i=0;i<NP;i++)
       {
           count++;
           CR=rand_normal(uCR,0.1,gen);
           F=rand_cauthy(uF,0.1,gen);
           /*
           for (int j=0;j<NP;j++)
           {
              P[j].fitness=calc_indi_fitness(P[j]);
           }
           */
           //cout<<"a"<<"    "<<calc_indi_fitness(P[0])<<endl;
           //cout<<"b"<<"    "<<P[0].fitness<<endl;
           sort(P.begin(),P.end());
           //sort(P.begin(),P.end(),cmp);
           //sort(P.begin(),P.begin()+NP,cmp);
           x_best=rand_int(0,(int)(p*NP)-1,gen);
           r1=rand_int(0,NP-1,gen);
           //while(r1==(xi=rand_int(1,NP,gen)));
           //xi=i+1;
           do
           {
              r1=rand_int(0,NP-1,gen);
           } while (r1==i);
 
           //cout<<"Asdfasdfsa"<<endl;
           do
           {
              r2=rand_int(0,NP+A.size()-1,gen);
           } while (r2==r1||r2==i);
           //cout<<count<<"  "<<x_best<<"  "<<xi<<"  "<<r1<<"  "<<r2<<endl;
           //if(r1==100)r1--;else if(r1==0)r1++;
           individual tmp;
           if(r2<NP)
              tmp=P[r2];
           else tmp=A[r2-NP];
           individual Vi,Ui;
 
           //xi--,x_best--,r1--;
           //Vi.fitness=P[xi].fitness+F*(P[x_best].fitness-P[xi].fitness+P[r1].fitness-tmp.fitness);
           for (int j=0;j<D;j++)
           {
              Vi.xreal[j]=P[i].xreal[j]+F*(P[x_best].xreal[j]-P[i].xreal[j]+P[r1].xreal[j]-tmp.xreal[j]);
           }
 
           //Vi=P[xi]+F*(P[x_best]-P[xi]+P[r1]-tmp);果然要重载
           int j_rand=rand_int(0,D-1,gen);
 
           for(int j=0;j<D;j++)
           {
              if(j==j_rand||rand_real(0,1,gen)<CR)
                  Ui.xreal[j]=Vi.xreal[j];
              else
                  Ui.xreal[j]=P[i].xreal[j];
           }
           double a=calc_indi_fitness(P[i]),b=calc_indi_fitness(Ui);
           if(a>b)
           if(calc_indi_fitness(P[i])>calc_indi_fitness(Ui))
           //{
           //  P[i]=P[i];
           //}else
           {
              //P[i]=Ui;
              for (int j=0;j<D;j++)
              {
                  P[i].xreal[j]=Ui.xreal[j];
              }
              A.push_back(P[i]);
              SF.push_back(F);
              SCR.push_back(CR);
           }
           P[i].fitness=a;
       }
       while (A.size()>NP)
       {
           int temp=rand_int(0,A.size()-1,gen);
           A.erase(find(A.begin(),A.end(),A[temp]));
       }
       uCR=(1.0-c)*uCR+c*Mean_A(SCR);
       uF=(1.0-c)*uF+c*Mean_L(SF);
       //cout<<calc_indi_fitness(P[0])<<endl;
       //cout<<P[0].fitness<<endl;
    }
    //sort(P.begin(),P.begin()+NP,cmp);
    sort(P.begin(),P.end());
    //cout<<calc_indi_fitness(P[0])<<endl;
    cout<<P[0].fitness<<endl;
    //int end_time=clock();
    //cout<<end_time-start_time<<endl;
    printf("%.6f\n",(double)clock()/CLOCKS_PER_SEC);
    system("pause");
    return 0;
}