1

这个问题可能很长，非常感谢您的耐心等待。核心问题是我使用 matlab 和 c++ 来实现优化算法，但它们为我提供了不同的结果（matlab 更好）。

我最近在研究一些进化算法，并对 PSO（粒子群优化）的一种变体感兴趣，它被称为竞争群优化器（2015 年出生）。这是论文链接http://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=6819057。

该算法的基本思想是首先在搜索空间中产生一些随机粒子并赋予它们随机速度。在每次迭代中，我们将它们随机配对，并让每对粒子比较它们的目标函数值。赢家（具有更好的客观价值观）保持现状，而输家则通过向赢家学习（走向赢家）来更新自己。

假设在迭代 t 时，粒子 i 和 j 被比较并且 i 更好。然后我们按照这些公式更新迭代 t+1 的粒子 j。如果粒子 j 超出搜索空间，我们只需将其拉回边界即可。R_1、R_2、R_3都是从[0, 1]均匀抽取的随机向量；操作“otimes”表示元素乘积；phi 是一个参数；x_bar 是 swarm 的中心。

例如，假设现在我要最小化一个 500 维的 Schwefel 函数（最小化最大绝对元素）并且我使用 250 个粒子，设置 phi=0.1，搜索空间为 500 维 [-100, 100]。当 C++ 卡在 85 到 90 之间时，Matlab 可以返回大约 35 的值。我无法弄清楚问题出在哪里。

让我在这里附上我的 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);
toc



function [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 number



double 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 的结果？非常感谢。

4

0 回答 0