参考:
http://blog.csdn.net/wangqiuyun/article/details/12515203
http://www.cnblogs.com/BreezeDust/p/3354769.html
第六章 群智能算法
《进化粒子群算法在TSP中的应用》
1.粒子群的全局版和局部版。全局版速度快,易陷入局部最优;局部版不易陷入局部最优,但速度慢。 具体可见 niuyongjie的专栏2.粒子群算法有多种方向。基本的粒子群算法可用来求函数最优值。 二进制粒子群算法 解决 组合优化问题(和遗传算法类似)。此外还有混合版等等。3.粒子群算法解决tsp问题。
tsp问题是一个组合优化问题,可以以类似遗传算法的方法解,参考 《进化粒子群算法在TSP中的应用》,我没有深入研究。还有一种速度用交换序来表示的方法,参考: http://blog.csdn.net/wangqiuyun/article/details/12515203, 第六章 群智能算法。我采用的是此法。
我们需要一个pbest来记录个体搜索到的最优解,用gbest来记录整个群体在一次迭代中搜索到的最优解。速度和粒子位置的更新公式如下:
v[i] = w * v[i] + c1 * rand() * (pbest[i] - present[i]) + c2 * rand() * (gbest - present[i])
present[i] = present[i] + v[i]
其中v[i]代表第i个粒子的速度,w代表惯性权值,c1和c2表示学习参数,rand()表示在0-1之间的随机数,pbest[i]代表第i个粒子搜索到的最优值,gbest代表整个集群搜索到的最优值,present[i]代表第i个粒子的当前位置。
本题在tsp问题上的算法描述:
Tsp问题的一个解为一个序列,可以表示为一个粒子,用一个序列的交换序列表示粒子的速度。
此时,速度和位置的更新公式为:
v[i] = w * v[i] + r1 * (pbest[i] - x[i]) + r2 * (gbest - x[i]);
x[i] = x[i] + v[i];
注:
w为权重值,设为0.729。
x[i]为粒子i的位置,pbeat[i]为粒子i的历史最佳位置,gbest[i]为全局历史最佳位置。
V[i]为粒子i的速度。
r1和r2是[0,1]的浮点数。
粒子群算法关键:
问题位置和速度的表示方法。
适度函数。
速度和位置的更新公式。(及参数的选取)
待加:
参考:百科,自话粒子群算法(超简单实例),等
标准粒子群算法解决最优值问题的实现?
参数的选择?
局部和全局的实现 ?(利用遗传算法或模拟退火解决陷入局部最优的问题)
一些参数选择:
粒子数: 可选30.一般取 20 – 40. 其实对于大部分的问题10个粒子已经足够可以取得好的结果, 不过对于比较难的问题或者特定类别的问题, 粒子数可以取到100 或 200
中止条件: 可选2000.最大循环数以及最小错误要求. 例如, 在上面的神经网络训练例子中, 最小错误可以设定为1个错误分类, 最大循环设定为2000, 这个中止条件由具体的问题确定.
学习因子: c1 和 c2 通常等于 2. 不过在文献中也有其他的取值. 但是一般 c1 等于 c2 并且范围在0和4之间
惯性权重:当Vmax很小时(对schaffer的f6函数,Vmax<=2),使用接近于1的惯性权重;当Vmax不是很小时(对schaffer的f6函数,Vmax>=3),使用权重w=0.8较好.如果没有Vmax的信息,使用0.8作为权重也是一种很好的选择.惯性权重w很小时偏重于发挥粒子群算法的局部搜索能力;惯性权重很大时将会偏重于发挥粒子群算法的全局搜索能力。
可选:0.729
全局PSO和局部PSO: 我们介绍了两种版本的粒子群优化算法: 全局版和局部版. 前者速度快不过有时会陷入局部最优. 后者收敛速度慢一点不过很难陷入局部最优. 在实际应用中, 可以先用全局PSO找到大致的结果,再用局部PSO进行搜索.
代码(不是很严格,可能有buge):
//#pragma warning (disable: 4786)
//#pragma comment (linker, "/STACK:16777216")
//HEAD
#include <cstdio>
#include <ctime>
#include <cstdlib>
#include <cstring>
#include <queue>
#include <string>
#include <set>
#include <stack>
#include <map>
#include <cmath>
#include <vector>
#include <iostream>
#include <algorithm>
#include <time.h>
#include <cstdlib>
using namespace std;
typedef long long LL;
const double MAX_VAL = (double)1e18;
const int MAX_GEN = 30;///最大迭代次数
const int MAX_SCALE = 3000;///最大种群规模
const int MAX_CITY = 20 + 2;///最大城市数
const double W_VAL = 0.729;///
struct SO{
int x, y;
SO(){}
SO(int x, int y): x(x), y(y){}
};
struct Point{
double x, y;
Point(){}
Point(int x, int y):x(x), y(y){};
void read()
{
scanf("%lf%lf", &x, &y);
}
};
inline int randomI(int x){ return rand()%x; }
inline double randomD(){ return (double)rand()/RAND_MAX; }
inline double getDist(Point a, Point b)
{
return sqrt((a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y));
}
struct PSO{
double w;
int scale;
int cityNum;
int nowGen;///当前代数
int maxGen;///迭代次数
int bestNum;
int bestGen;///最佳出现代数
double dist[MAX_CITY][MAX_CITY];
int oPop[MAX_SCALE][MAX_CITY];///粒子群
double fitness[MAX_SCALE];///种群适应度,表示种群中各个个体的适应度
vector<SO> listV[MAX_SCALE];/// 每科粒子的初始交换序列
int Pd[MAX_SCALE][MAX_CITY];///一颗粒子历代中出现最好的解,
double vPd[MAX_SCALE];///解的评价值
int Pgd[MAX_CITY];/// 整个粒子群经历过的的最好的解,每个粒子都能记住自己搜索到的最好解
double vPgd;/// 最好的解的评价值
PSO(){}
PSO(int s, int c, int mG, double ww, double d[MAX_CITY][MAX_CITY])
{
scale = s;
cityNum = c;
maxGen = mG;
w = ww;
for (int i = 0; i < cityNum; i++)
for (int j = 0; j < cityNum; j++)
dist[i][j] = d[i][j];
}
void copyArray(double a[], double b[], int n)
{
for (int i = 0; i < n; i++) a[i] = b[i];
}
void copyArray(int a[], int b[], int n)
{
for (int i = 0; i < n; i++) a[i] = b[i];
}
void init()
{
nowGen = 0;
for (int i = 0; i < scale; i++)
{
for (int j = 0; j < cityNum; )
{
int x = randomI(cityNum);
int r;
for (r = 0; r < j; r++)
{
if (x == oPop[i][r]) break;
}
if (r == j)
{
oPop[i][j] = x;
// cout << oPop[i][j] << ' ';
j++;
}
}
// cout << endl;
}
for (int i = 0; i < scale; i++)
{
// cout << i << " :" << endl;
int vn = randomI(cityNum) + 1;
for (int j = 0; j < vn; j++)
{
int x = randomI(cityNum);
int y = randomI(cityNum);
while (x == y) y = randomI(cityNum);
SO so(x, y);
listV[i].push_back(so);
// cout << so.x << "*" << so.y << ' ';
}
// cout <<endl;
}
getFitness();
for (int i = 0; i < scale; i++)
{
vPd[i] = fitness[i];
copyArray(Pd[i], oPop[i], cityNum);
}
bestNum = 0;
vPgd = fitness[0];
bestGen = 0;
for (int i = 0; i < scale; i++) if (vPgd > fitness[i])
{
vPgd = fitness[i];
bestNum = i;
}
copyArray(Pgd, oPop[bestNum], cityNum);
}
double getVal(int x)
{
double ret = 0;
for (int i = 0; i < cityNum; i++)
{
int xx = oPop[x][i % cityNum];
int yy = oPop[x][(i + 1) % cityNum];
ret += dist[xx][yy];
}
return ret;
}
void getFitness()
{
for (int i = 0; i < scale; i++)
fitness[i] = getVal(i);
}
void UpdateVal()
{
int j = 0;
double vj = fitness[0];
for (int i = 0; i < scale; i++)
{
if (vPd[i] > fitness[i])
{
vPd[i] = fitness[i];
copyArray(Pd[i], oPop[i], cityNum);///???
}
if (vj > fitness[i])
{
vj = fitness[i];
j = i;
}
}
if (vj < vPgd)
{
bestGen = nowGen;///
bestNum = j;///
vPgd = vj;
copyArray(Pgd, oPop[j], cityNum);
}
}
void changeTo(int a[], vector<SO> v)///
{
int vn = v.size();
for (int i = 0; i < vn; i++)
{
int x = v[i].x, y = v[i].y;
swap(a[x], a[y]);
}
}
vector<SO> minus(int a[], int b[])///
{
int c[MAX_CITY], d[MAX_CITY];
for (int i = 0; i < cityNum; i++) d[i] = b[i];
for (int i = 0; i < cityNum; i++) c[a[i]] = i;
vector<SO> v;
SO s;
for (int i = 0; i < cityNum; i++)
{
if (d[i] != a[i])
{
s.x = i, s.y = c[a[i]];
swap(d[s.x], d[s.y]);
v.push_back(s);
}
}
return v;
}
void addTo(vector<SO> &v, vector<SO> a, int vn)
{
for (int i = 0; i < vn; i++)
v.push_back(a[i]);
}
/// Vii=wVi+ra(Pid-Xid)+rb(Pgd-Xid)
void evolution()
{
for (int ig = 0; ig < maxGen; ig++)
{
nowGen = ig + 1;///nowGen
for (int is = 0; is < scale; is++)
{
if (is == bestNum) continue;
vector<SO> v;
v.clear();
int lvn = w * listV[is].size();
addTo(v, listV[is], lvn);
vector<SO> a = minus(Pd[is], oPop[is]);
int an = randomD() * a.size();
addTo(v, a, an);
vector<SO> b = minus(Pgd, oPop[is]);
int bn = randomD() * b.size();
addTo(v, b, bn);
listV[is] = v;
changeTo(oPop[is], listV[is]);
// cout << listV[is].size() << endl;
}
getFitness();
UpdateVal();
}
}
void solve()
{
init();
evolution();
printf("answer %lf:\n", vPgd);
printf("solution \n");
for (int i = 0; i < cityNum; i++)
{
if (i) printf(" ");
printf("%d\n", Pgd[i]);
}
cout << bestGen << ' ' << bestNum << endl;
puts("");
}
};
int cn;
vector<Point> pv;
double d[MAX_CITY][MAX_CITY];
void pre(vector<Point> pv)
{
int cn = pv.size();
for (int i = 0; i < cn; i++)
for (int j = i + 1; j < cn; j++)
{
d[i][j] = d[j][i] = getDist(pv[i], pv[j]);
}
}
int main ()
{
srand((unsigned long)time(0)); ///设置时间种子
while (cin >> cn)
{
pv.clear();
for (int i = 0; i < cn; i++)
{
Point p;
p.read();
pv.push_back(p);
}
pre(pv);
PSO solver(MAX_SCALE, cn, MAX_GEN, W_VAL, d);
cout << "***************************************************" << endl;
solver.solve();
cout << "***************************************************" << endl;
}
return 0;
}
// cout << nowGen << " ------------------------:" << endl;
// for (int i = 0; i < scale; i++)
// {
// cout << vPd[i] << ' ' << fitness[i] << endl;
// for (int j = 0; j < cityNum; j++)
// cout << Pd[i][j] << ' ';
// cout << endl;
// }
//
// for (int i = 0; i < scale; i++)
// {
// for (int j = 0; j < listV[i].size(); j++)
// cout << listV[i][j].x << "*" << listV[i][j].y << ' ';
// cout << endl;
// }
//
// cout << vPgd << ' ' << bestGen << ' ' << bestNum << ' ' << nowGen << endl;
// for (int i = 0; i < cityNum; i++) cout << Pgd[i] << ' ';
// cout << endl;
// cout << " ------------------------------------" << endl << endl;
/**
15
0 0
1 0
2 0
0 2
1 1
10 4
3 8
6 7
8 2
10 10
12 13
-10 1
2 5
1 100
12 12
*/