使用遗传算法计算方程的最大值
遗传算法介绍:
遗传算法使用计算机程序来模拟生物的进化过程,形象的来说,就是通过“繁殖”、“变异”,“竞争”以及自然选择来得到问题的最优解。
遗传算法的基本结构如下:
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满足:
(xmax- xmin) * 10t <= 2m - 1;
则由二进制编码(bin)到实数的转换如下:
x = xmin + decimal(bin) * (xmax- xmin ) / (2m - 1);
所以,在本例中,x可以用8位二进制串表示, y可以用5位二进制串表示。
确定评估函数:
由于此时就是计算这个方程的最大值,所以可以简单的取表达式f(x, y)来作为评估的函数。其函数值越大,表示评估的适应度越大。
确定选择策略:
选择的策略有很多。在本问题中,由于使用函数本身做评估函数,即评估的适应度越大,选择的概率应该越大,所以可以使用轮盘赌策略来进行选择,具体步骤如下:
计算种群中每个个体的适应度之和F = ∑eval(k);
计算每个个体的选择概率 Pk = eval(k) / F;
计算每个个体的累计概率 Qk = Qk-1 + Pk
产生一个0~1的随机数r, 如果Qk-1 <= r <= Qk , 则选择第k个个体。
确定遗传算子:
遗传算子主要指个体之间进行杂交和变异。
对于杂交,可以使用轮盘赌策略来选择部分父体进行信息的两两部分交换,交换的位置可以随机确定。
对于变异,同样可以使用轮盘赌策略来选择部分父体进行变异。随机的将表示这个个体的二进制位串的某个位求反。
终止条件:
由于使用函数本身做评估函数,故不可能根据评估的最大值来作为终止条件。这个时候可以选择迭代的次数作为终止的限制。
程序源码
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
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, double > > sol_set(COUNT); //最终选择个体的保存位置
vector< pair< bitset, 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 bitvec;
int j;
for(j = 0; j < BIT; j++)
...{
if(rand() % 2)
bitvec[j] = 1;
else
bitvec[j] = 0;
}
pair< bitset , double > sol(bitvec, 0.0);
sol_set[i] = sol;
}
}
//对种群中的个体进行评估
void evaluate()
...{
int i;
for(i = 0; i < COUNT; i++)
...{
bitset le(sol_set[i].first.to_string(), 0, LEFT_BIT);
bitset 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 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 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 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, double > sol;
bitset 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位)。计算的结果可能会接近这个值,也可能离这个值很远。多计算几次看看。