Contents

function gaSwarmExample
%GA_SWARM Demonstration of using GA to perform Swarm Optimization.
%
% Reference: http://en.wikipedia.org/wiki/Particle_swarm_optimization
%  Copyright (c) 2010, The MathWorks, Inc.
%  All rights reserved.

% Plot cost function isocontours
close all; clc
[XX,YY] = meshgrid(-5:.05:5); ZZ = zeros(size(XX));
for i = 1:size(XX,1)
    for j = 1:size(XX,2)
        ZZ(i,j) = rastriginsfcn([XX(i,j) YY(i,j)]);
    end
end
contourf(XX, YY, ZZ); colormap bone;
title('Rastrigin''s Function, Global Minimum at (0, 0)')
axis equal; axis tight;
hold on
plot(0, 0, 'r*');

% Turn off 'Elite' and 'Crossover' so that only mutations occur
% And all the mutations are simulations of a swarm
popSize = 100;    % Number of particles in the swarm
gaopts = gaoptimset( 'EliteCount', 0, 'CrossoverFraction', 0, ...
    'Generations', Inf, ...
    'PopulationSize', popSize, ...
    'InitialPopulation', 2*randn(popSize, 2) - 3, ...
    'MutationFcn', @mutationswarm, ...
    'SelectionFcn', @(~,nParents,~) (1:nParents) , ...
    'OutputFcn', @gaoutput );

% This is an unconstrained problem. Turn off constraints
A = []; b = []; Aeq = []; beq = [];
LB = []; UB = []; nonlcon = [];

% Call GA to perform Swarm Optimization through custom MutationFcn
ga(@rastriginsfcn, 2, A, b, Aeq, beq, LB, UB, nonlcon, gaopts)

end

MutationFcn to simulate Swarm Optimization

function mutationChildren = mutationswarm(~, ~, ~, ~, state, thisScore, ...
    thisPopulation, ~, ~)
persistent swarm_state      % Record state of the swarm
w = .95; C1 = 2; C2 = 2;    % Heuristic parameters can be adjusted

if state.Generation == 1    % Initialize state of the swarm
    swarm_state = [];
    swarm_state.X_lbest = zeros(size(thisPopulation));
    swarm_state.fitness_lbest = Inf(size(thisPopulation, 1), 1);
    swarm_state.X_gbest = zeros(1, size(thisPopulation, 2));
    swarm_state.fitness_gbest = Inf;
    swarm_state.V = zeros(size(thisPopulation));
end
mutationChildren = zeros(size(thisPopulation));

% For each individual to be mutated
for k = 1:size(thisPopulation, 1)

    % Update local best
    if swarm_state.fitness_lbest(k) > thisScore(k)
        swarm_state.fitness_lbest(k) = thisScore(k);
        swarm_state.X_lbest(k,:) = thisPopulation(k,:);
    end

end

% Update global best
[gbest, idx] = min(swarm_state.fitness_lbest);

if swarm_state.fitness_gbest > gbest
    swarm_state.fitness_gbest = gbest;
    swarm_state.X_gbest = thisPopulation(idx, :);
end

% Update velocity and position
for i = 1:size(thisPopulation, 1)        % For each individual
    for j = 1:size(thisPopulation, 2)    % Along each direction
        swarm_state.V(i,j) = w*swarm_state.V(i,j) ...
            + C1 * rand * (swarm_state.X_lbest(i,j) - thisPopulation(i,j)) ...
            + C2 * rand * (swarm_state.X_gbest(j) - thisPopulation(i,j)) ;
        mutationChildren(i,j) = thisPopulation(i,j) + 0.05*swarm_state.V(i,j);
    end
end
end

OutputFcn for GA

function [state,options,optchanged] = gaoutput(options,state,flag,~)
persistent l; optchanged = false;
if strcmpi(flag, 'init'); l = plot(0, 0, 'wo'); end
set(l, 'XData', state.Population(:,1), 'YData', state.Population(:,2));
pause(.05);
end