布谷鸟优化算法C++

#include <iostream>
#include <vector>
#include <cmath>
#include <random>
#include <time.h>
#include <fstream>
#define pi acos(-1)
//5只布谷鸟
constexpr int NestNum = 40;    
//pi值
//规定X,Y 的取值范围
constexpr double X_max = 5;
constexpr double X_min = 0;
constexpr double Y_max = 5;
constexpr double Y_min = 0;
//最大迭代次数
constexpr int MaxIterationTimes = 300;
//被宿主发现的概率
constexpr double Pa = 0.25;

//自变量结构体
struct Nest {
    double x;
    double y;
    double fitness;
};
void fitFunc(Nest& nest);
int findBetterNest(std::vector<Nest>&);
std::vector<Nest> levy(std::vector<Nest> OldNestPop, std::size_t bestNum);
std::vector<Nest> RandomAbandonPaNestPop(std::vector<Nest> OldNestPop);
//随机数引擎
static std::default_random_engine e(time(0));
static std::uniform_real_distribution<double> u(0, 1);
int main(void)
{
    bool flag_output = false;
    double Xold;
    double Xnew;
    double Yold;
    double Ynew;
    std::ofstream outfileX("D:\\cuckoo\\cuckooX.txt");
    std::ofstream outfileY("D:\\cuckoo\\cuckooY.txt");
    std::ofstream outfileZ("D:\\cuckoo\\cuckooZ.txt");
    //现在的鸟巢群
    std::vector<Nest> current_nestPop;
    //迭代次数
    int num = 0;
    //最优鸟巢
    int BestNestCurrent;    
    //初始化
    for (int i = 0; i < NestNum; ++i)
    {
        Nest nestinitial;
        nestinitial.x = (X_max - X_min) * u(e) + X_min;
        nestinitial.y = (Y_max - Y_min) * u(e) + Y_min;
        fitFunc(nestinitial);
        current_nestPop.push_back(nestinitial);
    }
    //for (auto i : nestPop)
    //{
    //    std::cout << i.fitness << std::endl;
    //}
    //寻找最优个体
    BestNestCurrent = findBetterNest(current_nestPop);
    outfileX << current_nestPop[BestNestCurrent].x << std::endl;
    outfileY << current_nestPop[BestNestCurrent].y << std::endl;
    outfileZ << current_nestPop[BestNestCurrent].fitness << std::endl;
    while (num < MaxIterationTimes)
    {
        //储存上次的最优解的X,Y
        Xold = current_nestPop[BestNestCurrent].x;
        Yold = current_nestPop[BestNestCurrent].y;
        //levy飞行--位置更新
        std::vector<Nest> NewNestPop = levy(current_nestPop, BestNestCurrent);
        //用适应值较好的鸟窝位置替换适应值较差的鸟窝位置
        for (decltype(NewNestPop.size()) i = 0; i < NewNestPop.size(); ++i)
        {
            if (i != BestNestCurrent && NewNestPop[i].fitness < current_nestPop[i].fitness)
            {
                current_nestPop[i] = NewNestPop[i];
            }
        }//此时得到更优的鸟窝位置
        //存安去险 保留鸟窝中被发现概率较小的鸟窝位置，并随机改变发现概率较大的鸟窝位置
        NewNestPop = RandomAbandonPaNestPop(current_nestPop);
        for (decltype(NewNestPop.size()) i = 0; i < NewNestPop.size(); ++i)
        {
            if (i != BestNestCurrent && NewNestPop[i].fitness < current_nestPop[i].fitness)
            {
                current_nestPop[i] = NewNestPop[i];
            }
        }//此时得到更优的鸟窝位置

        BestNestCurrent = findBetterNest(current_nestPop);//现在的最优鸟巢位置
        Xnew = current_nestPop[BestNestCurrent].x;
        Ynew = current_nestPop[BestNestCurrent].y;

        if (Xnew != Xold || Ynew != Yold)
        {
            outfileX << current_nestPop[BestNestCurrent].x << std::endl;
            outfileY << current_nestPop[BestNestCurrent].y << std::endl;
            
        }

        outfileZ << current_nestPop[BestNestCurrent].fitness << std::endl;
        /*std::cout << current_nestPop[BestNestCurrent].fitness << std::endl;
        std::cout << "(x,y)" << '(' << current_nestPop[BestNestCurrent].x << ',' << current_nestPop[BestNestCurrent].y << ')' << std::endl;*/
        //outfileX << current_nestPop[BestNestCurrent].x << std::endl;
        //outfileY << current_nestPop[BestNestCurrent].y << std::endl;
        //outfileZ << current_nestPop[BestNestCurrent].fitness << std::endl;

        num++;
    }

    std::cout << current_nestPop[BestNestCurrent].fitness << std::endl;
    return 0;
}

void fitFunc(Nest& nest)
{
    
    nest.fitness = -sin(nest.x) * pow(sin(nest.x * nest.x / pi), 20) - sin(nest.y) * pow(sin(2 * nest.y * nest.y / pi), 20);
    
    //nest.fitness = -(nest.x - 1) * (nest.x - 1) + 1;
}

int findBetterNest(std::vector<Nest>& nestPop)
{
    int BestNum = 0;
    for (decltype(nestPop.size()) i = 0; i < nestPop.size(); ++i)
    {
        if (nestPop[i].fitness < nestPop[BestNum].fitness)
        {
            BestNum = i;
        }
    }

    return BestNum;
}

std::vector<Nest> levy(std::vector<Nest> OldNestPop,std::size_t bestNum)
{
    double beta = 1.5;
    //    double alpha = 0.01 * R(e);//有的论文写
    double alpha = 0.4;
    double sigma_u = pow((tgamma(1 + beta) * sin(pi * beta / 2)) / (beta * tgamma((1 + beta) / 2) * pow(2, (beta - 1) / 2)), 1 / beta);
    double sigma_v = 1;

    static std::normal_distribution<double> R(0, sigma_u);
    static std::normal_distribution<double> R1(0, sigma_v);

    for (auto& i : OldNestPop) 
    {
        //前面的系数是保证最优鸟巢不会进行levy飞行
        double stepX = (i.x - OldNestPop[bestNum].x) * R(e) / (pow(abs(R1(e)), 1 / beta));
        double stepY = (i.x - OldNestPop[bestNum].x) * R(e) / (pow(abs(R1(e)), 1 / beta));
        //按范围更新X
        if (i.x + alpha * stepX > X_max)
        {
            i.x = X_max;
        }
        else if(i.x + alpha * stepX < X_min)
        {
            i.x = X_min;
        }
        else
        {
            i.x = i.x + alpha * stepX;
        }
        //按范围更新Y
        if (i.y + alpha * stepY > Y_max)
        {
            i.y = Y_max;
        }
        else if (i.y + alpha * stepY < Y_min)
        {
            i.y = Y_min;
        }
        else
        {
            i.y = i.y + alpha * stepY;
        }

        fitFunc(i);
    }

    return OldNestPop;
}

std::vector<Nest> RandomAbandonPaNestPop(std::vector<Nest> OldNestPop)
{
    double step_sizeX = 0;
    double step_sizeY = 0;
    static std::uniform_int_distribution<int> randomInt(0, OldNestPop.size() - 1);
    for(decltype(OldNestPop.size()) i = 0;i < OldNestPop.size();++i)
    {
        if (u(e) < Pa)//被宿主发现了，要重新寻找新巢
        {
            step_sizeX = u(e) * (OldNestPop[randomInt(e)].x - OldNestPop[randomInt(e)].x);
            step_sizeY = u(e) * (OldNestPop[randomInt(e)].y - OldNestPop[randomInt(e)].y);

            if (OldNestPop[i].x + step_sizeX > X_max)
            {
                OldNestPop[i].x = X_max;
            }
            else if(OldNestPop[i].x + step_sizeX < X_min)
            {
                OldNestPop[i].x = X_min;
            }
            else
            {
                OldNestPop[i].x += step_sizeX;
            }
            
            if (OldNestPop[i].y + step_sizeY > Y_max)
            {
                OldNestPop[i].y = Y_max;
            }
            else if (OldNestPop[i].y + step_sizeY < Y_min)
            {
                OldNestPop[i].y = Y_min;
            }
            else
            {
                OldNestPop[i].y += step_sizeY;
            }

            fitFunc(OldNestPop[i]);
        }
    }

    return OldNestPop;
}