source: liacs/nc/laser-pulse-shaping/pso.m@ 36

% Particle Swarm optimalisation
% BSDLicence
% Rick van der Zwet - 0433373 - <hvdzwet@liacs.nl>
% Modeled after http://en.wikipedia.org/wiki/Particle_swarm_optimization

% Dimention settings
parameters = 80;

% If global optimum does not change this many steps bail out 
iteration_break = 5;
max_iterations = 1000;
max_time = 5 * 60; % in sec

% Flock properties 
local_swarm_size = 10;
local_swarms = 50;

%% Particle properties
% Speed of walking around to a certain direction
wander = 0.4;

% 'Influence' of the envirionment with regards to solutions
% Trust the group global solution to be feasible
c_social  = 0.4;
% Trust the neighbor solution to be feasible
c_cognitive = 0.4;
% Trust the own best solution to be feasible
c_ego = 0.2;

% Variables used for plotting
fitness_history = [];
fitness_iterations = [];

% Initiate all particles
flock_p = rand(parameters,local_swarm_size,local_swarms) .* (2 * pi);
flock_v = zeros(size(flock_p));


% Global best placeholder
g_best = zeros(parameters,1);
g_fitness = 1;
% at (:,x) lives the neighbor best of local_swarm 'x'
n_best = zeros(parameters,local_swarms);
n_fitness = ones(parameters,local_swarms);
% at (:,p,x) leves the local best of particle 'p' in local_swarm 'x'
l_best = zeros(parameters,local_swarm_size,local_swarms);
l_fitness = ones(local_swarm_size, local_swarms);

idle_counter = 0;
start_time = time();

% Code not optimised for performance, but for readablility
for i = 1:max_iterations
    for s = 1:local_swarms
        fitness = SHGa(flock_p(:,:,s));
        % See if we got any better local optimum
        for p = 1:local_swarm_size
            if fitness(p) < l_fitness(p,s)
                l_fitness(p,s) = fitness(p);
                l_best(:,p,s) = flock_p(:,p,s);
            end
        end

        % See if we got any better neighbor optimum
        for p = 1:local_swarm_size
            if l_fitness(p,s) < n_fitness(s)
                n_fitness(s) = l_fitness(p,s);
                n_best(:,s) = l_best(:,p,s);
            end
        end
    end

    idle_counter = idle_counter + 1;

    % See wether we have a new global optimum
    for s = 1:local_swarms
        if n_fitness(s) < g_fitness
          g_fitness = n_fitness(s);
          g_best = n_best(:,s);
          idle_counter = 0;
        end
    end

    % Stop conditions
    if idle_counter == iteration_break 
        printf("Caught by idle_counter\n");
        break;
    end
    if time - start_time > max_time
        printf("Caught by max_time used \n");
        break;
    end
    

    printf("%i : %.15f\n", i, g_fitness);
    fitness_iterations = [fitness_iterations, i];
    fitness_history = [fitness_history, g_fitness];
    

    % Update particles to new value
    r_cognitive = rand();
    r_social = rand();
    r_ego = rand();
    for s = 1:local_swarms
        for p = 1:local_swarm_size
            flock_v(:,p,s) = flock_v(:,p,s) * wander;
            flock_v(:,p,s) = flock_v(:,p,s) + (g_best - flock_p(:,p,s)) * (c_cognitive * r_cognitive);
            flock_v(:,p,s) = flock_v(:,p,s) + (n_best(:,s) - flock_p(:,p,s)) * (c_social * r_social);
            flock_v(:,p,s) = flock_v(:,p,s) + (l_best(:,p,s) - flock_p(:,p,s)) * (c_ego * r_ego);
            flock_p(:,p,s) = flock_p(:,p,s) + flock_v(:,p,s);
        end
    end
end


plot(fitness_iterations,fitness_history);
title(sprintf('Particle Swarm Optimalisation on Laser-Pulse shaping problem'));
ylabel('fitness');
xlabel('iterations');
grid on;
legend(sprintf('Parameters %i',parameters));
print("pso-fitness.eps","-depsc2");