使用遗传算法计算方程的最大值

• 遗传算法介绍：

       遗传算法使用计算机程序来模拟生物的进化过程，形象的来说，就是通过“繁殖”、“变异”，“竞争”以及自然选择来得到问题的最优解。

      遗传算法的基本结构如下：

 

t  =   0 ;
initialize(p(t));      // 找出原始解集
evaluate(p(t));    // 评估解集中是否有最优解
while (not termination condition)
do
    pp(t)  =  parent_select(p(t));     // 选择父体
    pp(t)  =  crossover(pp(t));         // 父体之间互换信息
    p(t + 1 )  =  mutate(pp(t));           // 变异，引入多样性
    evaluate(p(t + 1 ));                   // 评估是否有最优解
    t  =  t + 1 ;
end  while

return  the best solution

 

• 具体实例

 

• •  问题描述：

                           求解问题：max f(x, y) = 21.5 + x * sin(4x * PI) + y * sin(20y * PI); 

                           其中 -3.0 <= x <= 12.1, 4.1 <= y <=5.8; 要求解精确到小数点之后1位。

• • 个体表示：

如果我们要用遗传算法解决问题，第一步要解决的是如何表示每一个个体。一般情况下，我们使用一个二进制的整数来表示一个个体，整数的位数与所要表示的精度有关，如：

如果用m表示x(精确到小数点后t位)所需要的二进制位串的长度，则m满足：

(x_max- x_min) * 10^t <= 2^m - 1;

 则由二进制编码(bin)到实数的转换如下：

x ＝ x_min + decimal(bin) * (x_max- x_min ) / (2^m - 1);

所以，在本例中，x可以用8位二进制串表示， y可以用5位二进制串表示。

• 确定评估函数： 

由于此时就是计算这个方程的最大值，所以可以简单的取表达式f(x, y)来作为评估的函数。其函数值越大，表示评估的适应度越大。

•  确定选择策略：

选择的策略有很多。在本问题中，由于使用函数本身做评估函数，即评估的适应度越大，选择的概率应该越大，所以可以使用轮盘赌策略来进行选择，具体步骤如下：

1. 计算种群中每个个体的适应度之和F ＝  ∑eval(k)；
2. 计算每个个体的选择概率 P _k = eval(k) / F；
3. 计算每个个体的累计概率 Q _k ＝ Q _k－1 + P _k
4. 产生一个0～1的随机数r， 如果Q _k－1 <= r  <=  Q _k ， 则选择第k个个体。

• 确定遗传算子：

遗传算子主要指个体之间进行杂交和变异。

对于杂交，可以使用轮盘赌策略来选择部分父体进行信息的两两部分交换，交换的位置可以随机确定。

对于变异，同样可以使用轮盘赌策略来选择部分父体进行变异。随机的将表示这个个体的二进制位串的某个位求反。

• 终止条件：

由于使用函数本身做评估函数，故不可能根据评估的最大值来作为终止条件。这个时候可以选择迭代的次数作为终止的限制。

• 程序源码

 

#include  < iostream >
#include  < vector >
#include  < utility >
#include  < bitset >
#include  < algorithm >
#include  < string >
#include  < cmath >
#include  < cstdlib >
#include  < ctime >
#include  < cassert >
using   namespace  std;

const   double  PI  =   3.1415926535 ;

const   double  pc  =   0.75 ;      // 杂交概率
const   double  pm  =   0.015 ;   // 变异概率

const   double  X_LOW  =   - 3.0 ;
const   double  X_UP  =   12.1 ;
const   double  Y_LOW  =   4.1 ;
const   double  Y_UP  =   5.8 ;

const   short  LEFT_BIT  =   8 ;
const   short  RIGH_BIT  =   5 ;
const   short  BIT  =  LEFT_BIT  +  RIGH_BIT;

const   short  GEN  =   1000 ;    // 迭代次数
const   short  COUNT  =   20 ;   // 种群规模
vector <  pair <  bitset < BIT > ,  double   >   >  sol_set(COUNT);   // 最终选择个体的保存位置
vector <  pair <  bitset < BIT > ,  double   >   >  temp(COUNT);

void  initialize();
void  evaluate();
void  parent_select();
void  crossover();
void  mutate();

// 目标函数
double  func( double  x,  double  y)
... {
    return 21.5 + x * sin(4 * x * PI) + y * sin(20 * PI * y);
}
// 使用随机数初始化种群
void  initialize()
... {
    srand((unsigned short) time(NULL));
    int i;
    for(i = 0; i < COUNT; i++)
    ...{
        bitset<BIT> bitvec;
        int j;
        for(j = 0; j < BIT; j++)
        ...{
            if(rand() % 2)
                bitvec[j] = 1;
            else
                bitvec[j] = 0;
        }
        
        pair< bitset<BIT> , double > sol(bitvec, 0.0);
        sol_set[i] = sol;
    }
}

// 对种群中的个体进行评估
void  evaluate()
... {
    int i;
    for(i = 0; i < COUNT; i++)
    ...{
        bitset<LEFT_BIT> le(sol_set[i].first.to_string(), 0, LEFT_BIT);
        bitset<RIGH_BIT> ri(sol_set[i].first.to_string(), LEFT_BIT, RIGH_BIT);
        
        double x = X_LOW + le.to_ulong() * (X_UP - X_LOW) / ((1 << LEFT_BIT) - 1);
        assert(x >= X_LOW && x <= X_UP);

        double y = Y_LOW + ri.to_ulong() * (Y_UP - Y_LOW) / ((1 << RIGH_BIT) - 1);
        assert(y >= Y_LOW && y <= Y_UP);

        sol_set[i].second =  func(x, y);
    }
}
// 使用轮盘赌策略选择个体
void  parent_select()
... {
    double sum = 0.0, p = 0.0, q = 0.0;
    int i;
    for(i = 0 ; i < COUNT; i++)
    ...{
        sum += sol_set[i].second;
    }

    vector<double> pq(COUNT);
    for(i = 0 ; i < COUNT; i++)
    ...{
        p = sol_set[i].second / sum;
        q += p;
        pq[i] = q;
    }

    srand((unsigned short) time(NULL));
    for(i = 0 ; i < COUNT; i++)
    ...{
        double r = rand() * 1.0 / RAND_MAX;

        int j;
        for (j = 0; r > pq[j] ; j++);
        
        temp[i] = sol_set[j];
    }
    sol_set.swap(temp);
}
// 按照轮盘赌策略杂交
void  crossover()
... {
    srand((unsigned short) time(NULL));
    vector<int> index;
    int i ;
    for(i = 0; i < COUNT; i++)
    ...{
        double r = rand() * 1.0 / RAND_MAX;
        if( r < pc)
        ...{
            index.push_back(i);
        }
    }

    int m, n, j;
    while(index.size() >= 2)
    ...{
        m = rand() % index.size();
        index.erase(index.begin() + m);

        n = rand() % index.size();
        index.erase(index.begin() + n);

        i = rand() % BIT;
        
        for(j = i; j < BIT; j++)
        ...{
            if(sol_set[m].first[j] != sol_set[n].first[j])
            ...{
                sol_set[m].first.flip(j);
                sol_set[n].first.flip(j);
            }
        }
    }    
}
// 按照轮盘赌策略变异
void  mutate()
... {
    srand((unsigned short) time(NULL));
    vector<int> index;
    int i ;
    for(i = 0; i < COUNT; i++)
    ...{
        double r = rand() * 1.0 / RAND_MAX;
        if( r < pm)
        ...{
            index.push_back(i);
        }
    }

    pair< bitset<BIT>, double > sol;
    bitset<BIT> a;

    while(index.size() > 0)
    ...{
        sol = sol_set.at(index.at(index.size() - 1));

        i = rand() % BIT;

        a = sol.first;
        a.flip(i);

        index.pop_back();
    }
}

int  main( int  argc,  char   ** argv)
... {
    int i;
    initialize();
    evaluate();

    for(i = 0 ; i < GEN; i++)
    ...{
        parent_select();
        crossover();
        mutate();
        evaluate();
    }

    for(i = 0; i < COUNT; i++)
    ...{
        cout << sol_set.at(i).second << endl;
    }
    return 0;
}

 

程序比较简单，基本上都是自解释的。

由于遗传算法经过大量的随即选择过程，故最终的结果不一定正确。在这个题目中，按照微积分的方法算出最大值为38.5左右（精确到小数点后1位）。计算的结果可能会接近这个值，也可能离这个值很远。多计算几次看看。