Matlab и C ++ дают разные результаты при оптимизации функции Швефеля по алгоритму, аналогичному PSO

Этот вопрос может быть длинным, и я очень ценю ваше терпение. Основная проблема заключается в том, что я использовал matlab и c ++ для реализации алгоритма оптимизации, но они дали мне разные результаты (лучше matlab).

Недавно я изучал некоторые эволюционные алгоритмы и заинтересовался одним вариантом PSO (Particle Swarm Optimization), который называется Competitive Swarm Optimizer (родился в 2015 году). Это бумажная ссылка http://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber = 6819057.

Основная идея этого алгоритма состоит в том, чтобы сначала генерировать случайные частицы в пространстве поиска и назначать им случайные скорости. На каждой итерации мы случайным образом соединяем их и позволяем каждой паре частиц сравнивать значения их целевых функций. Победители (с лучшими объективными значениями) сохраняют статус-кво, в то время как проигравшие обновляют себя, учась у победителей (переходя к победителям).

Предположим, что на итерации t частицы i и j сравниваются, и i лучше. Затем мы обновляем частицу j для итерации t + 1, следуя этим формулам. Если частица j находится вне пространства поиска, мы просто возвращаем ее обратно к границе. R_1, R_2, R_3 — все случайные векторы, равномерно взятые из [0, 1]; операция «otimes» означает поэлементное произведение; фи является параметром; x_bar это центр роя.

Например, предположим, что теперь я хочу минимизировать 500-дневную функцию Швефеля (минимизировать максимальный абсолютный элемент), и я использую 250 частиц, установите phi = 0,1, пространство поиска равно 500-д [-100, 100]. Matlab мог вернуть мне что-то около 35, в то время как C ++ застрял в 85 до 90. Я не могу понять, в чем проблема.

Позвольте мне присоединить мой код Matlab и C ++ здесь.

Sch = @(x)max(abs(x))
lb = -100 * ones(1, 500);
ub = 100 * ones(1, 500);
swarmsize = 250;
phi = 0.1;
maxiter = 10000;
tic
cso(Sch, lb, ub, swarmsize, phi, maxiter);
tocfunction [minf, minx] = cso(obj_fun, lb, ub, swarmsize, phi, maxiter)
assert(length(lb) == length(ub), 'Not equal length of bounds');
if all(ub - lb <= 0) > 0
error('Error. \n Upper bound must be greater than lower bound.')
end
vhigh = abs(ub - lb);
vlow = -vhigh;
S = swarmsize;
D = length(ub);
x = rand(S, D);
x = bsxfun(@plus, lb, bsxfun(@times, ub-lb, x)); % randomly initalize all particles
v = zeros([S D]); % set initial velocities to 0
iter = 0;
pairnum_1 = floor(S / 2);
losers = 1:S;
fx = arrayfun(@(K) obj_fun(x(K, :)), 1:S);
randperm_index = randperm(S);
while iter <= maxiter
fx(losers) = arrayfun(@(K) obj_fun(x(K, :)), losers);
swarm_center = mean(x); % calculate center all particles
randperm_index = randperm(S); % randomly permuate all particle indexes
rpairs = [randperm_index(1:pairnum_1); randperm_index(S-pairnum_1+1:S)]'; % random pair
cmask= (fx(rpairs(:, 1)) > fx(rpairs(:, 2)))';
losers = bsxfun(@times, cmask, rpairs(:, 1)) + bsxfun(@times, ~cmask, rpairs(:, 2)); % losers who with larger values
winners = bsxfun(@times, ~cmask, rpairs(:, 1)) + bsxfun(@times, cmask, rpairs(:, 2)); % winners who with smaller values
R1 = rand(pairnum_1, D);
R2 = rand(pairnum_1, D);
R3 = rand(pairnum_1, D);
v(losers, :) = bsxfun(@times, R1, v(losers, :)) + bsxfun(@times, R2, x(winners, :) - x(losers, :)) + phi * bsxfun(@times, R3, bsxfun(@minus, swarm_center, x(losers, :)));
x(losers, :) = x(losers, :) + v(losers, :);
maskl = bsxfun(@lt, x(losers, :), lb);
masku = bsxfun(@gt, x(losers, :), ub);
mask = bsxfun(@lt, x(losers, :), lb) | bsxfun(@gt, x(losers, :), ub);
x(losers, :) = bsxfun(@times, ~mask, x(losers, :)) + bsxfun(@times, lb, maskl) + bsxfun(@times, ub, masku);
iter = iter + 1;
fprintf('Iter: %d\n', iter);
fprintf('Best fitness: %e\n', min(fx));
end
fprintf('Best fitness: %e\n', min(fx));
[minf, min_index] = min(fx);
minx = x(min_index, :);
end

(Я не писал функцию C ++.)

#include <cstring>
#include <iostream>
#include <cmath>
#include <algorithm>
#include <ctime>
#include <iomanip>
#include <time.h>
#include <math.h>
#include <stdlib.h>
#include <stdio.h>

#define rand_01 ((double) rand() / RAND_MAX)  // generate 0~1 random numbers
#define PI 3.14159265359

const int numofdims = 500;  // problem dimension
const int numofparticles = 250;  // number of particles
const int halfswarm = numofparticles / 2;
const double phi = 0.1;
const int maxiter = 10000; // iteration numberdouble Sch(double X[], int d); // max(abs(x_i))using namespace std;

int main(){
clock_t t1,t2;
t1=clock();

srand(time(0)); // random seed

double** X = new double*[numofparticles];  // X for storing all particles
for(int i=0; i<numofparticles; i++)
X[i] = new double[numofdims];

double** V = new double*[numofparticles]; // V for storing velocities
for(int i=0; i<numofparticles; i++)
V[i] = new double[numofdims];

double Xmin[numofdims] = {0}; // lower bounds
double Xmax[numofdims] = {0}; // upper bounds

double* fitnesses = new double[numofparticles];  // objective function values

for(int j=0; j<numofdims; j++)
{
Xmin[j] = -100;
Xmax[j] = 100;
}

for(int i=0; i<numofparticles; i++)
{
for(int j=0; j<numofdims; j++)
{
X[i][j] = Xmin[j] + rand_01 * (Xmax[j] - Xmin[j]);  // initialize X
V[i][j] = 0; // initialize V
}
}

for(int i=0; i<numofparticles; i++)
{
fitnesses[i] = Sch(X[i], numofdims); //
}

double minfit = fitnesses[0]; // temporary minimal value
int minidx = 0; // temporary index of minimal value

int* idxofparticles = new int[numofparticles];

for(int i=0; i<numofparticles; i++)
idxofparticles[i] = i;

double* Xmean = new double[numofdims];

int* losers = new int[halfswarm]; // for saving losers indexes

for(int iter=0; iter<maxiter; iter++)
{

random_shuffle(idxofparticles, idxofparticles+numofparticles);

for(int j=0; j<numofdims; j++)
{
for(int i=0; i<numofparticles; i++)
{
Xmean[j] += X[i][j];
}
Xmean[j] = (double) Xmean[j] / numofparticles;  // calculate swarm center
}

for(int i = 0; i < halfswarm; i++)
{
// indexes are now random
// compare 1st to (halfswarm+1)th, 2nd to (halfswarm+2)th, ...
if(fitnesses[idxofparticles[i]] < fitnesses[idxofparticles[i+halfswarm]])
{

losers[i] = idxofparticles[i+halfswarm];
for(int j = 0; j < numofdims; j++)
{
V[idxofparticles[i+halfswarm]][j] = rand_01 * V[idxofparticles[i+halfswarm]][j] + rand_01 * (X[idxofparticles[i]][j] - X[idxofparticles[i+halfswarm]][j]) + rand_01 * phi * (Xmean[j] - X[idxofparticles[i+halfswarm]][j]);
X[idxofparticles[i+halfswarm]][j] = min(max((X[idxofparticles[i+halfswarm]][j] + V[idxofparticles[i+halfswarm]][j]), Xmin[j]), Xmax[j]);
}
}
else
{
losers[i] = idxofparticles[i];
for(int j = 0; j < numofdims; j++)
{
V[idxofparticles[i]][j] = rand_01 * V[idxofparticles[i]][j] + rand_01 * (X[idxofparticles[i+halfswarm]][j] - X[idxofparticles[i]][j]) + rand_01 * phi * (Xmean[j] - X[idxofparticles[i]][j]);
X[idxofparticles[i]][j] = min(max((X[idxofparticles[i]][j] + V[idxofparticles[i]][j]), Xmin[j]), Xmax[j]);
}
}

}
// recalculate particles' values
for(int i=0; i<numofparticles; i++)
{
fitnesses[i] = Sch(X[i], numofdims);
if(fitnesses[i] < minfit)
{
minfit = fitnesses[i]; // update minimum
minidx = i; // update index
}
}if(iter % 1000 == 0)
{
cout << scientific << endl;
cout << minfit << endl;

}
}
cout << scientific << endl;
cout << minfit << endl;

t2=clock();

delete [] X;
delete [] V;
delete [] fitnesses;
delete [] idxofparticles;
delete [] Xmean;
delete [] losers;

float diff ((float)t2-(float)t1);
float seconds = diff / CLOCKS_PER_SEC;
cout << "runtime: " << seconds << "s" <<endl;

return 0;
}double Sch(double X[], int d)
{
double result=abs(X[0]);
for(int j=0; j<d; j++)
{
if(abs(X[j]) > result)
result = abs(X[j]);
}
return result;
}

Итак, наконец, почему мой код на c ++ не может воспроизвести результат matlab? Большое спасибо.