多目标 NSGAII

这其实是上个礼拜就完成的了，但由于上个礼拜没有开会

这周三才开的会，然后确认了算法的正确性了，今天有时间就来记录下

这个也是基于IEEE的论文写的

花的时间比较多，而且，也是我写过的最大的一个算法了

先用我自己的语言来描述下整体的算法吧：

一开始当然是初始化种群P了（在这里可以同时初始化一个种群Q）

先把P和Q归并到R

然后构造边界集F，也就是快速非支配的排序（这是个重点以及难点）

构造完了新的边界集之后就要产生新的种群P了

这一步过程中，要对每一层进行拥挤距离的计算（为了保持特种的多样性）

拥挤距离的计算也是个难点+重点

对于需要计算的最后一层，还要进行一层偏序关系的排序来进行选取

P构造好了，就要通过P生成Q了

这里，首先要采用二元锦标赛选择产生NP/2个个体（NP是种群大小）

然后通过杂交变异生成新的个体（随机数小于0.9进行SBX，即杂交，否则变异）

杂交是两个个体产生两个个体的，变异是一个个体产生一个个体

当Q的数量达到了NP之后，第一轮迭代结束

整个算法流程大概是这样吧，最后我跑出的结果还行

收敛和分布都还可以，还要等用matlab绘图之后看最后的结果了

下面就直接贴代码了

这次为了自己能记清（太大了，思路总不清晰），我大部分地方都加了注释了

#include<cstdio>
#include<cstring>
#include<cmath>
#include<cmath>
#include<ctime>
#include<iostream>
#include<cstdlib>
#include<algorithm>
using namespace std;

const int NP=200;
const int INF=~0U>>2;
const double upper=1.0;
const double lower=0.0;
const int D=30;
const int G=2000;
const int problem=2;
const double cc=20.0;
const double mm=20.0;

struct individual{
    double xreal[D+10];
    double fitness[problem+2];
    int num_bedominated;
    int set_dominated[2*NP+10];
    int num_dominated;
    int rank;
    double distance;
    friend bool operator < (const individual&a,const individual&b){
        return a.distance>b.distance;
    }
};
individual P[NP+10];
individual Q[NP+10];
individual R[2*NP+10];
struct Front{
    individual ind[2*NP+10];
    int num;
}F[2*NP+10];
void output_pop(individual P[]){
    for(int i=0;i<NP;i++){

        if(P[i].xreal[0]>1.0){
            cout<<i<<"     ";
        for(int j=0;j<5;j++)
            cout<<P[i].xreal[j]<<"   ";
        cout<<endl;
        }
    }
}
void init_xreal(individual &P){
    for(int i=0;i<D;i++)
        P.xreal[i]=rand()/(RAND_MAX+0.0);
}
void calc_fitness_ZDT1(individual &P){
    P.fitness[0]=P.xreal[0];
    double tmp=0.0;
    for(int i=1;i<D;i++)
        tmp+=P.xreal[i];
    double g_x=1.0+9*tmp/(D-1);
    P.fitness[1]=g_x*(1.0-sqrt(P.xreal[0]/g_x));
    /*if(P.fitness>1.0){
    ;
    }*/
    //printf("--%f  %f\n",P.fitness[0],P.fitness[1]);
}
void init_population(individual pop[],int N){
    for(int i=0;i<N;i++){
        pop[i].distance=0.0;
        init_xreal(pop[i]);
        calc_fitness_ZDT1(pop[i]);
    }
}
void mergePQ(individual R[],individual P[],individual Q[]){
    for(int i=0;i<NP;i++)
        R[i]=P[i];
    for(int i=NP;i<2*NP;i++)
        R[i]=Q[i-NP];
}
bool is_dominated(individual pop[],int p,int q){
    if(pop[p].fitness[0]<=pop[q].fitness[0]&&
        pop[p].fitness[1]<=pop[q].fitness[1]){
            if(pop[p].fitness[0]<pop[q].fitness[0]||
                pop[p].fitness[1]<pop[q].fitness[1])
                return true;
    }
    return false;
}
int fast_nondominated_sort(individual pop[],int N){
    int F_num=0;
    int number[2*NP+10]={0};
    for(int p=0;p<N;p++){
        int num_dominated=0;
        int num_bedominated=0;
        for(int q=0;q<N;q++){
            if(p==q)continue;
            if(is_dominated(pop,p,q)){
                pop[p].set_dominated[num_dominated++]=q;
            }else
                if(is_dominated(pop,q,p)){
                    num_bedominated++;
                }
        }
        number[p]=num_bedominated;

        pop[p].num_dominated=num_dominated;

        pop[p].num_bedominated=num_bedominated;
        if(num_bedominated==0){
            pop[p].rank=0;
            F[0].ind[F_num++]=pop[p];
        }
    }
    F[0].num=F_num;
    int m=0;
    while(F[m].num!=0){
        F_num=0;
        for(int p=0;p<F[m].num;p++){
            //printf("p:  %d\n",p);
            individual &ind=F[m].ind[p];
            //printf("num_domi:  %d\n",ind.num_dominated);
            for(int q=0;q<ind.num_dominated;q++){
                int nq=pop[ ind.set_dominated[q] ].num_bedominated;
                nq--;
                int id=ind.set_dominated[q];
                number[id]--;
                nq=number[id];
                if(nq==0){
                    pop[ ind.set_dominated[q] ].rank=m+1;
                    F[m+1].ind[F_num++]=pop[ ind.set_dominated[q] ];
                }
            }
        }
        F[++m].num=F_num;
    }
    return m;
}
bool cmp0(individual a,individual b){
    return a.fitness[0]<b.fitness[0];
}
bool cmp1(individual a,individual b){
    return a.fitness[1]<b.fitness[1];
}
void crowding_distance_assignment(individual pop[],int N){
    for(int i=0;i<N;i++)pop[i].distance=0;
    for(int pro=0;pro<problem;pro++){
        switch(pro){
        case 0:sort(pop,pop+N,cmp0);break;
        case 1:sort(pop,pop+N,cmp1);break;
        }
        pop[0].distance=pop[N-1].distance=INF+0.0;
        for(int i=1;i<N-1;i++)
            pop[i].distance += ( pop[i+1].fitness[pro] - pop[i-1].fitness[pro] ) / ( pop[N-1].fitness[pro] - pop[0].fitness[pro] );
    }
}
bool cmp(individual a,individual b){
    if(a.rank<b.rank)
        return true;
    else
        if(a.rank==b.rank && a.distance>b.distance)
            return true;
    return false;
}
void BitwiseMutation_for_BinaryCoded(const individual P[],individual tmp[]){
    int index,a,b;
    bool flag[NP+10]={0};
    for(int i=0;i<NP/2;i++){
        while(true){
            a=rand()%NP;
            if(flag[a])continue;
            break;
        }
        while(true){
            b=rand()%NP;
            if(b==a)continue;
            if(flag[b])continue;
            break;
        }
        if(cmp(P[a],P[b]))
            index=a;
        else
            index=b;
        flag[index]=true;
        tmp[i]=P[index];
    }
}
double judge(double x){
    if(x<0.0)
        x=0.0;
    else
        if(x>1.0)
            x=1.0;
    return x;
}
void SBX_(individual Q[],individual tmp[],int&j){
    double u,beta;
    int a,b;
    a=rand()%(NP/2);
    while(true){
        b=rand()%(NP/2);
        if(a==b)continue;
        break;
    }
    individual inda=tmp[a];
    individual indb=tmp[b];
    for(int i=0;i<D;i++){
        //while(true){
            u=rand()/(RAND_MAX+1.0);
            //if(u==0.0)continue;
            //break;
        //}
        if(u<0.5)
            beta=pow(2*u,1/(cc+1));
        else
            beta=1/pow(2*(1-u),1/(cc+1));
        double y1=0.5*((1-beta)*inda.xreal[i]+(1+beta)*indb.xreal[i]);
        double y2=0.5*((1+beta)*inda.xreal[i]+(1-beta)*indb.xreal[i]);
        //cout<<"y1:   "<<y1<<"   y2:   "<<y2<<endl;
        Q[j].xreal[i]=judge(y1);
        Q[j+1].xreal[i]=judge(y2);
    }
    calc_fitness_ZDT1(Q[j]);
    calc_fitness_ZDT1(Q[j+1]);
    //j+=2;
}
void polynomial_mutation(individual Q[],individual tmp[],int&j){
    int a=rand()%(NP/2);
    double rk,beta;
    individual ind=tmp[a];
    double rate=1.0/D;

    for(int i=0;i<D;i++){
        double random=rand()/(RAND_MAX+1.0);
        if(random>=rate)
        {
            Q[j].xreal[i]=tmp[a].xreal[i];
            continue;
        }
        rk=rand()/(RAND_MAX+1.0);

        if(rk<0.5)
            beta=pow(2*rk,1/(mm+1))-1.0;
        else
            beta=1.0-pow(2*(1.0-rk),1/(mm+1));
        double y=ind.xreal[i]+(upper-lower)*beta;
        Q[j].xreal[i]=judge(y);
    }
    calc_fitness_ZDT1(Q[j]);
    //j++;
}

void mamk_new_pop(const individual P[],individual Q[]){
    individual tmp[NP];
    BitwiseMutation_for_BinaryCoded(P,tmp);
    for(int i=0;i<NP;i++){
        double r=rand()/(RAND_MAX+0.0);
        //cout<<r<<endl;
        if(r<0.9){
            SBX_(Q,tmp,i);
        }else
            polynomial_mutation(Q,tmp,i);
    }
}


int main(){
    freopen("nsga2.out","w",stdout);
    srand((unsigned int)time(NULL));
    init_population(P,NP);
    init_population(Q,NP);

    for(int gen=0;gen<G;gen++){
        //printf("gen:  %d\n",gen);
        mergePQ(R,P,Q);

        int m=fast_nondominated_sort(R,2*NP);

        int num_P=0;
        int i=0;
        while(num_P+F[i].num<=NP){
            crowding_distance_assignment(F[i].ind,F[i].num);
            //printf("%d\n",F[i].num);

            for(int j=0;j<F[i].num;j++)
                P[num_P++]=F[i].ind[j];

            //printf("num_P:  %d\n",num_P);
            i++;
        }
        //printf("gen1:  %d\n",gen);

        //crowding_distance_assignment(F[i].ind,F[i].num);
        if(num_P<NP){
            crowding_distance_assignment(F[i].ind,F[i].num);
            //cout<<F[i].num<<endl;
            //sort(F[i].ind,F[i].ind+F[i].num,cmp);
            sort(F[i].ind,F[i].ind+F[i].num);
            int tmp=num_P;
            for(int j=0;j<NP-tmp;j++)
                P[num_P++]=F[i].ind[j];
        }
        //cout<<num_P<<endl;
        //output_pop(P);
        mamk_new_pop(P,Q);
        //output_pop(P);
    }

    for(int i=0;i<NP;i++)printf("%f %f\n",P[i].fitness[0],P[i].fitness[1]);
    /*puts("");
    for(int i=0;i<NP;i++)calc_fitness_ZDT1(P[i]);
    for(int i=0;i<NP;i++)printf("%f %f\n",P[i].fitness[0],P[i].fitness[1]);*/

    for(int i=0;i<NP;i++){
        cout<<i<<"     ";
        for(int j=0;j<D;j++)
            cout<<P[i].xreal[j]<<"   ";
        cout<<endl;
    }
    printf("time:  %.6f\n",(double)clock()/CLOCKS_PER_SEC);
    //system("pause");
    return 0;
}