CMAGABC算法
`// ABC.cpp : Defines the entry point for the console application.
//自动测30个函数，输出加上成功率
include
include
include
include <stdlib.h>
include <math.h>
include <time.h>
include
include
using namespace std;
const int NP = 50; //种群规模
const int FoodNumber = NP / 2; //蜜源数量=引领蜂数量（每个蜜源同一时间内只有一只蜜蜂采蜜）跟随蜂数量也是NP/2
const int D = 30; //维数30
const int limit = 200; //经过t次迭代到达limit没找到更好的蜜源，引领蜂转化为侦查蜂
const int maxCycle = 100* D;
const int runtime = 50; //最长运行时间
const int maxFES= 10000 * D; //最大迭代次数
//const double PI = 3.14159265358979323846; //pai值
const double Mt=0.1;
const int Gen = maxFES / 100; // Output Gen. 50 (D = 30) 75 (D = 60) 125 (D = 100)
int best; //当前最好的解
double Foods[FoodNumber][D];
double f[FoodNumber];
double trial[FoodNumber];
double prob[FoodNumber];
double solution[D];
double ObjValSol0;
double ObjValSol;
double FitnessSol;
double GlobalMin;
double GlobalMins[runtime];
double lb;
double ub;
double sigma;
double weights[FoodNumber]; //中间值
double weights2[FoodNumber]; //系数wi
int mu; //系数μ
double m[D]; //参数m(g+1)
double m_old[D]; //保存上一代mg
double cc;
double cs;
double mucov;
double mueff;
double ccov;
double damps;
double muti[D][D]; //矩阵相乘的中间值
int FEs, FEsmin, FEsmean, acount, SucNum;
double mean, Std, SucRate, Accept = 0.0; // 1.0e-8; //
typedef double (*FunctionCallback)(double sol[D]);
//#include "dfs.h"
include "cec2014.h" // D = 10 40 * 2500 = 100000 100 * 1000 D = 30 40 * 7500 = 300000 100 * 3000
FunctionCallback function;//另一种，一起调用三十个
typedef struct { //排序数组返回索引
int index;
double value;
}sort_ar;
bool compare(sort_ar a,sort_ar b){
return a.value<b.value; //比较函数，用于sort函数中，升序排列
}
void initRandom(int index) //初始化 侦查蜂xi的初始化
{
int j;
double r;
r = (double)rand()/(double)(RAND_MAX) ; //第一个随机数有可能一样（初始化r）
for(j=0;j<D;j++)
{
r = (double)rand()/(double)(RAND_MAX) ;
Foods[index][j]=r*(ub-lb)+lb; //给每个Food[i][j]赋值 侦查蜂公式
}
f[index]=function( Foods[index] ); FEs++; //迭代次数加一
trial[index]=0; //迭代次数回到初始状态0
}
void initial() //给每个蜜源赋值
{
int i;
for(i=0;i<FoodNumber;i++)
{
initRandom(i);
}
GlobalMin = f[0];
}
void CalculateProbabilities() //计算适应值
{
int i;
double sum, tempfun;
double fitness[FoodNumber];
for (i=0;i<FoodNumber;i++)
{
tempfun = f[i]; //fi
if( tempfun >= 0 )
fitness[i] = 1.0 / ( 1.0 + tempfun ); //fi>=0的适应值 1/（1+fi）
else
fitness[i] = 1.0 - tempfun; //fi<=0的适应值 1+abs（fi）
}
sum=fitness[0];
for (i=1;i<FoodNumber;i++)
{
sum=sum+fitness[i];
}
for (i=0;i<FoodNumber;i++)
{
prob[i]=fitness[i]/sum; //pi=fitnessi/fitness（1到SN）的和
}
}
int Roulette( int k = 0 ) //轮盘赌
{
int i, j;
double r, m;
r = (double)rand()/(double)(RAND_MAX) ; //r的初始化，在0,1内随机
if ( k == 0 ){
r = (double)rand()/(double)(RAND_MAX) ; //r在0,1内随机
k = (int)(r*FoodNumber) % FoodNumber;
}
r = (double)rand()/(double)(RAND_MAX) ; //r在0,1内随机
m = 0;
for(i=0; i<FoodNumber; i++)
{
j = ( i + k ) % FoodNumber;
m = m + prob[j]; //m为适应值的累积
if( r < m ) //r<m 在i周围产生新蜜源
{
break;
}
}
if( j >= FoodNumber ){
r = (double)rand()/(double)(RAND_MAX) ;
j = (int)(r*FoodNumber) % FoodNumber;
}
return j; //轮盘赌的结果返回j
}
void MemorizeBestSource() //记住当前最优解
{
int i;
double temp;
temp = f[0];
for(i=1;i<FoodNumber;i++)
{
if (f[i] < temp){
temp=f[i];
best=i;
}
}
if(GlobalMin > temp)
GlobalMin = temp;
}
double pc[D][1] = {0}; //协方差矩阵C的进化路径
double ps[D][1] = {0}; //sigma步长的进化路径
double Bmartix[D][D] = {0};
double Dmartix[D][D] = {0};
double Cmartix[D][D] = {0}; //协方差矩阵
double D_1martix[D][D] = {0}; //D的逆矩阵
double B_tmartix[D][D] = {0}; //B的转置
double c[D][D] = {0}; //C的负二分之一次方
double chiN;
double arz[D][FoodNumber] = {0};
void CMAInit(){
int i,j;
double sum = 0.0;
sigma = 0.5;
mu = floor(FoodNumber / 2); //精英集大小
for(i = 0;i < mu;i++){
weights[i] = log(mu + 1) - log(i+1);
sum += weights[i];
}
for(i = 0;i<mu;i++){
weights2[i] = weights[i] / sum; //求wi
}
sum = 0.0;
for(i = 0;i<mu;i++){
sum += pow(weights2[i],2);
}
mueff = 1 / sum; //mueff
cc = 4.0 /(D+4.0);
cs = (mueff + 2.0) / (D + mueff + 3.0);
mucov = mueff;
ccov = (1.0 / mucov) * 2.0 / pow((D + 1.4),2) + (1.0 - 1.0/mucov) * ((2.0 * mueff - 1.0) / (pow(D + 2.0 , 2) + 2.0 * mueff ));
damps = 1.0 + 2.0 * max(0.0,sqrt((mueff-1) / (D+1.0)) - 1.0) + cs;
for(i = 0;i < D;i++){ //生成单位矩阵
Bmartix[i][i] = 1;
Dmartix[i][i] = 1;
Cmartix[i][i] = 1;
m[i] = ((double)rand() / (double)(RAND_MAX));
ps[i][0] = 0;
pc[i][0] = 0;
}
chiN = pow(D,0.5) * (1.0 - 1.0 / (4.0 * D) + 1.0 / (21.0 * pow(D,2)));
}
bool mutiply(double *a,double *b,int dim,double * mut){
int i,j,k;
for(i = 0;i < dim;i++){
for(j = 0;j < dim;j++){
mut[i*dim+j] = 0.0;
for(k=0; k < dim;k++){
mut[i*dim+j] += a[i*dim+k] * b[k*dim + j];
}
}
}
return true;
}
//求向量的2范数
double norm(double arc[D][1]){
double ans = 0;
for(int i=0;i<D;i++){
ans += pow(arc[i][0],2);
}
ans = sqrt(ans);
return ans;
}
//求对角矩阵的逆矩阵
bool DiagonalMatrixInverse(double arc[D][D],int n,double invarc[D][D]){
int i;
for(i = 0;i<n;i++){
invarc[i][i] = 1.0/arc[i][i];
}
return true;
}
//矩阵转置
bool transposition(double arc[D][D],int n,double trarc[D][D]){
int i,j;
for(i = 0;i<n;i++){
for(j = 0;j<n;j++){
trarc[i][j] = arc[j][i];
}
}
return true;
}
/雅克比迭代法特征分解，matrix为目标矩阵，dim矩阵行列数，precision精度，max最大迭代次数，matrix为特征值矩阵
/
bool Jacobi(double matrix, int dim, double eigenvectors, double precision, int maxtime)
{
for (int i = 0; i < dim; i++) { //单位矩阵
eigenvectors[idim + i] = 1.0f;
for (int j = 0; j < dim; j++) {
if (i != j)
eigenvectors[idim + j] = 0.0f;
}
}
int nCount = 0; //current iteration
while (1) {
//找到矩阵上三角部分绝对值最大的元素的行列数
double dbMax = matrix[1];
int nRow = 0;
int nCol = 1;
for (int i = 0; i < dim; i++) { //row
for (int j = 0; j < dim; j++) { //column
double d = fabs(matrix[i*dim + j]);
if ((i != j) && (d > dbMax)) {
dbMax = d;
nRow = i;
nCol = j;
}
}
}
if (dbMax < precision) //precision check
break;
if (nCount > maxtime)
break; //iterations check
nCount++;
double dbApp = matrix[nRow*dim + nRow];
double dbApq = matrix[nRow*dim + nCol];
double dbAqq = matrix[nCol*dim + nCol];
//compute rotate angle
double dbAngle = 0.5*atan2(-2.0 * dbApq, dbAqq - dbApp);
double dbSinTheta = sin(dbAngle);
double dbCosTheta = cos(dbAngle);
double dbSin2Theta = sin(2.0 * dbAngle);
double dbCos2Theta = cos(2.0 * dbAngle);
matrix[nRow*dim + nRow] = dbApp*dbCosTheta*dbCosTheta +
dbAqq*dbSinTheta*dbSinTheta + 2 * dbApq*dbCosTheta*dbSinTheta;
matrix[nCol*dim + nCol] = dbApp*dbSinTheta*dbSinTheta +
dbAqq*dbCosTheta*dbCosTheta - 2 * dbApq*dbCosTheta*dbSinTheta;
matrix[nRow*dim + nCol] = 0.5*(dbAqq - dbApp)*dbSin2Theta + dbApq*dbCos2Theta;
matrix[nCol*dim + nRow] = matrix[nRow*dim + nCol];
for (int i = 0; i < dim; i++) {
if ((i != nCol) && (i != nRow)) {
int u = i*dim + nRow; //p
int w = i*dim + nCol; //q
dbMax = matrix[u];
matrix[u] = matrix[w] * dbSinTheta + dbMax*dbCosTheta;
matrix[w] = matrix[w] * dbCosTheta - dbMax*dbSinTheta;
}
}
for (int j = 0; j < dim; j++) {
if ((j != nCol) && (j != nRow)) {
int u = nRow*dim + j; //p
int w = nCol*dim + j; //q
dbMax = matrix[u];
matrix[u] = matrix[w] * dbSinTheta + dbMax*dbCosTheta;
matrix[w] = matrix[w] * dbCosTheta - dbMax*dbSinTheta;
}
}
//compute eigenvector
for (int i = 0; i < dim; i++) {
int u = i*dim + nRow; //p
int w = i*dim + nCol; //q
dbMax = eigenvectors[u];
eigenvectors[u] = eigenvectors[w] * dbSinTheta + dbMax*dbCosTheta;
eigenvectors[w] = eigenvectors[w] * dbCosTheta - dbMax*dbSinTheta;
}
}
return true;
}
double diff[D][D]; //m(g+1)-m(g)
sort_ar sar[FoodNumber];
void CMAgenerate(){
int i,j,k;
double diff_x[D][mu]; //xg+1-mg
double diff_y[mu][D]; //上面的转置
double hsig;
for(i = 0;i < FoodNumber;i++){ //排序专用结构体
sar[i].index = i;
sar[i].value = f[i];
}
sort(sar,sar+FoodNumber,compare); //algorithm库里的排序函数
for(i = 0;i < mu;i++){
for(j = 0;j < D;j++){
m[j] += weights2[i] * Foods[sar[i].index][j];
}
}
for(j = 0;j < D;j++){
diff[j][0] = m[j] - m_old[j];
}
bool re = DiagonalMatrixInverse(Dmartix,D,D_1martix); //求D的逆和B的转置
re = transposition(Bmartix,D,B_tmartix);
re = mutiply(&Bmartix[0][0],&D_1martix[0][0],D, &muti[0][0]); //五行算Cg^-1/2 * m(g+1)-mg
re = mutiply(&muti[0][0],&B_tmartix[0][0],D,&c[0][0]);
re = mutiply(&c[0][0],&diff[0][0],D,&muti[0][0]);
for(j = 0;j < D;j++){
ps[j][0] = (1 - cs) * ps[j][0] + sqrt(cs * (2-cs) * mueff) * muti[j][0] / sigma;
}
hsig = norm(ps)/sqrt(1 - pow(1-cs,(2 * FEs / FoodNumber))/ chiN < 1.4+2/(D+1));
for(j = 0;j < D;j++){
pc[j][0] = (1 - cc) * pc[j][0] + hsig * sqrt(cc * (2-cc) * mueff) * diff[j][0] / sigma;
}
for(i = 0;i < mu;i++){
for(j = 0;j < D;j++){
diff_x[j][i] = (Foods[sar[i].index][j] - m[j]) / sigma;
diff_y[i][j] = diff_x[j][i];
}
}
//sigma = sigma * exp((cs/damps)*(norm(ps)/chiN - 1));
for(i = 0;i < D;i++){
for(j = 0;j < D;j++){
muti[i][j] = 0.0;
for(k=0; k < mu;k++){
muti[i][j] += diff_x[i][k] * diff_y[k][j];
}
}
}
//更新协方差矩阵
for(i = 0;i < D;i++){
for(j = 0;j < D;j++){
c[i][j] = (1.0-ccov) * c[i][j] + ccov * ((1.0/mucov) * pc[i][0] * pc[j][0] + (1.0 - hsig) * cc * (2.0 - cc) * c[i][j]) + ccov * (1.0 - 1.0/mucov) * weights2[i] * muti[i][j] ;//pc[i][0]*pc[j][0]就是伪代码中的pc*pc'
Dmartix[i][j] = c[i][j];
}
}
re = Jacobi(&Dmartix[0][0], D ,&Bmartix[0][0] , 1e-4,1000);
for(i = 0;i < D;i++){
for(j = 0;j < D;j++){
if(i != j) Dmartix[i][j] = 0.0;
}
Dmartix[i][i] = sqrt(abs(Dmartix[i][i]));
m_old[i] = m[i];
m[i] = 0;
}
}
double gaussrand(double E, double V) //高斯分布随机数
{
static double V1,V2,S;
static int phase = 0;
double X;
if(phase == 0){
do {
double U1 = (double)rand()/RAND_MAX;
double U2 = (double)rand()/RAND_MAX;
V1 = 2 * U1 - 1;
V2 = 2 * U2 - 1;
S = V1 * V1 + V2 * V2;
}while(S >= 1 || S == 0);
X = V1 * sqrt(-2 * log(S) / S);
} else {
X = V2 * sqrt(-2 * log(S) / S);
}
phase = 1 - phase;
X = X * V + E;
return X;
}
void SendEmployedBeesStd() //雇佣蜂阶段
{
int i,j;
int neighbour, dim2change;
double r;
for(i=0;i<FoodNumber;i++)
{
r = (double)rand()/(double)(RAND_MAX) ;
dim2change=(int)(r * D) % D;
r = (double)rand()/(double)(RAND_MAX) ;
neighbour =(int)(r * FoodNumber) % FoodNumber;
while(neighbour == i)
{
r = (double)rand()/(double)(RAND_MAX) ;
neighbour =(int)(r*FoodNumber) % FoodNumber;
}
for(j=0;j<D;j++)
{
solution[j]=Foods[i][j];
}
ObjValSol0=f[i];
r = 2.0 * (double)rand()/(double)RAND_MAX - 1.0;
/*for(j=0;j<D;j++){
solution[j]=Foods[i][j]+(Foods[i][j]-Foods[neighbour][j])*r ; //vij
if ( solution[j] < lb)
solution[j] = lb;
if ( solution[j] > ub)
solution[j] = ub;
}*/
solution[dim2change]=Foods[i][dim2change]+(Foods[i][dim2change]-Foods[neighbour][dim2change])*r ; //vij
while ( solution[dim2change] < lb)
solution[dim2change] = lb;
while (solution[dim2change] > ub)
solution[dim2change] = ub;
ObjValSol=function(solution); FEs++;
if (ObjValSol<ObjValSol0) //如果vi的值大于xi
{
trial[i]=0;
/*for(j=0;j<D;j++){
Foods[i][j]=solution[j]; //那么xi=vi
}*/
Foods[i][dim2change]=solution[dim2change]; //那么xi=vi
f[i]=ObjValSol; SucNum++;//
if( ObjValSol <= Accept ){
FEsmin = FEs;
break;
}
} else {
trial[i]=trial[i]+1; //迭代次数加一
}
}
}
void SendOnlookerBeesStd() //跟随蜂阶段
{
int i,j,t;
int neighbour, dim2change;
double r,sum;
bool re;
for(t=0;t<FoodNumber;t++)
{
sum = 0;
i = Roulette( ); //i=轮盘赌的结果j
r = (double)rand()/(double)(RAND_MAX) ;
dim2change=(int)(r*D) % D;
r = (double)rand()/(double)(RAND_MAX) ;
neighbour=(int)(r * FoodNumber) % FoodNumber;
while(neighbour == i)
{
r = (double)rand()/(double)(RAND_MAX) ;
neighbour=(int)(r*FoodNumber) % FoodNumber;
}
re = mutiply(&Bmartix[0][0],&Dmartix[0][0],D,&muti[0][0]);
arz[dim2change][i] = gaussrand(0.0,1.0);
for(j=0;j<D;j++)
{
solution[j]=Foods[i][j];
}
ObjValSol0=f[i];
r = 2.0 * (double)rand()/(double)RAND_MAX - 1.0;
//solution[dim2change]=Foods[i][dim2change]+(Foods[i][dim2change]-Foods[neighbour][dim2change])*r ;
for(j = 0;j < D;j++){
sum += muti[dim2change][j] * arz[dim2change][i];
}
solution[dim2change] = m_old[dim2change] + sigma * sum; //x取样
while ( solution[dim2change]<lb)
solution[dim2change] =lb;
while (solution[dim2change]>ub)
solution[dim2change] = ub;
/* while ( solution[dim2change] < lb )
solution[dim2change] += ub - lb;
while (solution[dim2change] > ub)
solution[dim2change] += lb - ub;
*/
ObjValSol=function(solution); FEs++;
if (ObjValSol<ObjValSol0)
{
trial[i]=0;
Foods[i][dim2change]=solution[dim2change];
f[i]=ObjValSol; SucNum++;//
if( ObjValSol <= Accept ){
FEsmin = FEs;
break;
}
} else {
trial[i]=trial[i]+1; //迭代次数加一
}
}
}
void SendScoutBees() //侦查蜂阶段
{
int i;
for(i=0;i<FoodNumber;i++)
{
if( trial[i] > limit ) //i的迭代次数大于limit的话，引领蜂转化为侦查蜂，随机产生新的xi
{
initRandom(i);
}
}
}
void ABC(int Fun)
{
//FunctionCallback functionArr[2] = { &f1, &f2 }; //, &f2, &f3, &f4, &f5, &f6, &f7, &f8, &f9, &f10, &f11, &f12, &f13, &f14, &f15,&f16, &f17, &f18, &f19, &f20, &f21, &f22, &f23, &f24, &f25, &f26, &f27, &f28, &f29, &f30
int iter, run, count, index;
struct Result {
int FEs;
double fit;
};
Result Red[100];
fstream ffit;
int i;
mean=0; Std=0; index = 1; SucRate = 0;
FEsmean = 0;
count = 0;
cout << "Limit = " << limit << endl;
srand(time(NULL));
function( solution );
for(run=0; run<runtime; run++)
{
initial();
CMAInit();
FEs = 0; FEsmin = 0; SucNum = 0;
for (iter=0;iter<maxCycle;iter++)
{
SendEmployedBeesStd( );
if( FEsmin != 0 )
break;
MemorizeBestSource();
CalculateProbabilities();
SendOnlookerBeesStd(); // Onlookbees Algorithm
if( FEsmin != 0 )
break;
MemorizeBestSource();
SendScoutBees();
CMAgenerate();
if( run == 3 && FEs / Gen == index ){
Red[index-1].FEs = index * Gen;
Red[index-1].fit = GlobalMin;
index++;
}
}
MemorizeBestSource();
if( run == 3 && FEs / Gen == index ){
Red[index-1].FEs = index * Gen;
Red[index-1].fit = GlobalMin;
index++;
}
cout.width(2);
cout << run + 1 << ".run: " << setiosflags(ios::scientific) << setprecision(2) << " GlobalMin = " << GlobalMin << " Num: " << best + 1
<< " FEsMin: " << FEsmin << " FEs: " << FEs << " Success rate: " << double( SucNum ) / double( FEs ) << endl;
GlobalMins[run]=GlobalMin;
mean=mean+GlobalMin;
SucRate = SucRate + double( SucNum ) / double( FEs );
if( GlobalMin <= Accept ){ // *** result ***
FEsmean = FEsmean + FEsmin;
acount++;
}
}
/*
ffit.open( "fitfile.dat", ios::out|ios::binary );
ffit.write( (char *)Red, sizeof( Result ) * maxCycle /Gen );
ffit.close( );
cout << endl;
for( i = 0; i < maxCycle /Gen; i++ ){
cout.width(2);
cout << Red[i].cycle << " " << Red[i].fit << " ";
cout << endl;
}
cout << endl;
*/
/*
// *** Evolutionary process data collection EPDC ***
char Filename[30] = "CABC_Test_Fun_";
char buf[10];
itoa( Fun, buf, 10 );
strcat( Filename, buf );
strcat( Filename, ".txt" );
ffit.open( Filename, ios::out);
ffit << " *** CABC Algorithm Test Result ***" << endl << endl;
ffit << "Function " << Fun << endl;
ffit << setiosflags(ios::scientific) << setprecision(2);
ffit << "FEs " << "GlobalMin " << endl;
for( i = 0; i < index - 1; i++ ){
ffit << Red[i].FEs << " " << Red[i].fit << endl;
}
ffit << endl;
ffit.close( );
cout << endl;
// *** Evolutionary process data collection END ***
*/
if( count != 0 )
FEsmean = FEsmean / acount;
mean = mean / runtime;
SucRate = SucRate / runtime;
for(run = 0; run < runtime; run++) {
Std = Std + pow(GlobalMins[run]-mean, 2);
}
Std = sqrt(Std/runtime);
cout << "Means of " << runtime << " runs: " << setiosflags(ios::scientific) << setprecision(2) << mean << endl;
cout << "Std of " << runtime << " runs: " << setiosflags(ios::scientific) << setprecision(2) << Std << endl;
cout << "SucRate = " << SucRate << endl;
cout << "Means of FEs = " << FEsmean << endl;
cout << "The number of accept (0.00) is "<< acount << endl;
}
int main()
{
FunctionCallback functionArr[30] = { &f1, &f2, &f3, &f4, &f5, &f6, &f7, &f8, &f9, &f10, &f11, &f12, &f13, &f14, &f15,
&f16, &f17, &f18, &f19, &f20, &f21, &f22, &f23, &f24, &f25, &f26, &f27, &f28, &f29, &f30 };
fstream ffit;
ffit.open( "ABC30.txt", ios::out );
ffit << " *** ABC Algorithm Test Result ***" << endl << endl;
ffit << "Parameter setting: " << "NP = " << NP << " D = " << D << " MaxFES = " << maxFES << " runtime = " << runtime << endl << endl;
int i;
for( i = 0; i < 30; i++){
function = functionArr[i];
ABC(i + 1);
fileflag = 0;
ffit << "Function :"<< i+1 << " Means of " << runtime << " runs: " << setiosflags(ios::scientific) << setprecision(2) << mean;
ffit << " Std of " << runtime << " runs: " << Std << setiosflags(ios::scientific) << setprecision(4) << " SucRate = " << SucRate << endl;
ffit << setiosflags(ios::scientific) << setprecision(2) << "Mean( Std ): "<< mean << "("<< Std << ")";
ffit << " Means of FEs = " << FEsmean << " The number of accept(0.00) is "<< acount << endl << endl;
}
ffit.close();
return 0;
}
`