粒子群优化算法(Particle Swarm Optimization, PSO)代码template (matlab)

% set the target function
objectiveFunction = @yourtargetfunction;

options.NumParticles            = yourparticlenum;
options.MinNeighborhoodFraction = 0.4;
options.seed                    = 11222;
options.inertia_range           = [0.1 0.4];
options.cSelf                   = 1;
options.cSocial                 = 1;
options.StallIterLimit          = 5;
options.MaxFunEvals             = 3000;
options.MaxIter                 = 25;
options.ObjectiveLimit          = 0.0001;
options.TolFun                  = 1.0E-4;

inertia_range = options.inertia_range;
cSelf         = options.cSelf;
cSocial       = options.cSocial;
min_inertia   = min(inertia_range);
max_inertia   = max(inertia_range);
maxInertia    = max_inertia;
minInertia    = min_inertia;

numParticles  = options.NumParticles;
seed          = options.seed;

minNeighborhoodSize = max(2,floor(numParticles*0.4));
rng(seed);

% initiate the parameter vector
paras_init = your initial parameter vector; 
numVariables = length(paras_init);

% set upper bound and lower bound
for i = 1:numVariables
    lb(i) = use your knowledge to determine it;
    ub(i) = use your knowledge to determine it;
end

lb_mat        = repmat(reshape(lb,1,numVariables),numParticles,1);
ub_mat        = repmat(reshape(ub,1,numVariables),numParticles,1);

% The initial position, which will be used as the P position later
p_position = rand(numParticles,numVariables);
p_position = p_position.*lb_mat + (1-p_position).*ub_mat;

% Enforce the bounds
p_position = max(lb_mat,p_position);
p_position = min(ub_mat,p_position);

values      = zeros(numParticles,1);

for i = 1:numParticles
    fprintf(1,'-----------------------\n');
    fprintf('Particle %d is Initiating Started\n', i);
    fprintf(1,'%-30s\n', 'The Position is')
    p_position(i,:)
    values(i) = objectiveFunction(p_position(i,:));
    fprintf(1,'------------------------\n');
end
    
% --------------------------------------------------
% The initial state;
% --------------------------------------------------
state.Fvals       = values;
state.Positions   = p_position;

state.IndividualBestFvals     = state.Fvals;
state.IndividualBestPositions = state.Positions;
state.FunEval                 = numParticles;

[bestFvals,bestPos]  = min(state.Fvals);
fprintf(1,'%-30s:%-10g\n','Best Fvals',bestFvals);
bestFvalsWindow = nan(options.StallIterLimit, 1);

vmax = ub - lb;

v    = repmat(-vmax,numParticles,1) + ...
       repmat(2*vmax,numParticles,1) .* rand(numParticles,numVariables);
   
state.Velocities = v;

inertia = max(inertia_range);

adaptiveInertiaCounter = 0;

if all(inertia_range >= 0)
    adaptiveInertia = max_inertia;
elseif all(inertia_range <= 0)
    adaptiveInertia = min_inertia;
end

adaptiveNeighborhoodSize = minNeighborhoodSize;
state.Iteration = 0;

exitFlag = [];
while isempty(exitFlag)
    
    fprintf(1,'--------------------------------------------------\n');
    fprintf('Iteration %d Started\n', state.Iteration);
    state.Iteration = state.Iteration + 1;
    
    neighborIndex = zeros(numParticles, adaptiveNeighborhoodSize);
    neighborIndex(:, 1) = 1:numParticles; % First neighbor is self
    
    for i = 1:numParticles
        % Determine random neighbors that exclude the particle itself,
        % which is (numParticles-1) particles
        neighbors = randperm(numParticles-1, adaptiveNeighborhoodSize-1);
        % Add 1 to indicies that are >= current particle index
        iShift = neighbors >= i;
        neighbors(iShift) = neighbors(iShift) + 1;
        neighborIndex(i,2:end) = neighbors;
    end
    % Identify the best neighbor
    [~, bestRowIndex] = min(state.IndividualBestFvals(neighborIndex), [], 2);
    bestLinearIndex = (bestRowIndex.'-1).*numParticles + (1:numParticles);
    bestNeighborIndex = neighborIndex(bestLinearIndex);
    
    % Randomness in velocities
    randSelf = rand(numParticles, numVariables);
    randSocial = rand(numParticles, numVariables);
    
    newVelocities = adaptiveInertia*state.Velocities + ...
    cSelf*randSelf.*(state.IndividualBestPositions-state.Positions) + ...
    cSocial*randSocial.*(state.IndividualBestPositions(bestNeighborIndex, :)-state.Positions);
    tfValid = all(isfinite(newVelocities), 2);
    
    % Check for illegal values
    state.Velocities(tfValid,:) = newVelocities(tfValid,:);
    
    % Update the positions
    newPopulation = state.Positions + state.Velocities;
    tfInvalid = ~isfinite(newPopulation);
    newPopulation(tfInvalid) = state.Positions(tfInvalid);
    
    % Enforce bounds, setting the corresponding velocity component to
    % zero if a particle encounters a lower/upper bound
    tfInvalid = newPopulation < lb_mat;
    if any(tfInvalid(:))
        newPopulation(tfInvalid) = lb_mat(tfInvalid);
        state.Velocities(tfInvalid) = 0;
    end
    tfInvalid = newPopulation > ub_mat;
    if any(tfInvalid(:))
        newPopulation(tfInvalid) = ub_mat(tfInvalid);
        state.Velocities(tfInvalid) = 0;
    end
    state.Positions = newPopulation;
    
    values      = zeros(numParticles,1);
    tic;
    for i = 1:numParticles
        fprintf(1,'%-30s:%-10d;%-30s:%-10d\n','Current Iteration:',state.Iteration,'Current Particle',i);
        values(i) = objectiveFunction(p_position(i,:));
    end
    time = toc;
    fprintf(1,'%-30s:%-30g\n','Time for eval',time);
    state.FunEval = state.FunEval  + numParticles;

    state.Fvals = values;
    % Remember the best fvals and positions
    tfImproved = state.Fvals < state.IndividualBestFvals;
    state.IndividualBestFvals(tfImproved) = state.Fvals(tfImproved);
    state.IndividualBestPositions(tfImproved, :) = state.Positions(tfImproved, :);
    bestFvalsWindow(1+mod(state.Iteration-1,options.StallIterLimit)) = min(state.IndividualBestFvals);
    
    % Keep track of improvement in bestFvals and update the adaptive
    % parameters according to the approach described in S. Iadevaia et
    % al. Cancer Res 2010;70:6704-6714 and M. Liu, D. Shin, and H. I.
    % Kang. International Conference on Information, Communications and
    % Signal Processing 2009:1-5.
    newBest = min(state.IndividualBestFvals);
    if isfinite(newBest) && newBest < bestFvals
        bestFvals = newBest;
        state.LastImprovement = state.Iteration;
        adaptiveInertiaCounter = max(0, adaptiveInertiaCounter-1);
        adaptiveNeighborhoodSize = minNeighborhoodSize;
    else
        adaptiveInertiaCounter = adaptiveInertiaCounter+1;
        adaptiveNeighborhoodSize = min(numParticles, adaptiveNeighborhoodSize+minNeighborhoodSize);
    end
    
    % Update the inertia coefficient, enforcing limits (Since inertia
    % can be negative, enforcing both upper *and* lower bounds after
    % multiplying.)
    if adaptiveInertiaCounter < 2
        adaptiveInertia = max(minInertia, min(maxInertia, 2*adaptiveInertia));
    elseif adaptiveInertiaCounter > 5
        adaptiveInertia = max(minInertia, min(maxInertia, 0.5*adaptiveInertia));
    end
    fprintf(1,'%-30s:%-10g\n','Best Fvals',bestFvals);
    [exitFlag,reasonToStop] = stopPso(options,state,bestFvalsWindow);
    fprintf(1,'--------------------------------------------------\n');
end

if exitFlag == 1 || exitFlag == -1 
        % Find and return the best solution
        [fval,bestFvals] = min(state.IndividualBestFvals);
        x = state.IndividualBestPositions(bestFvals,:);

        % Update output structure
        output  = state;
        
        positionBest = state.IndividualBestPositions;
        fvalBest = state.IndividualBestFvals;
        
end

function [exitFlag,reasonToStop] = stopPso(options,state,bestFvalsWindow)
    iteration = state.Iteration;

    iterationIndex = 1+mod(iteration-1,options.StallIterLimit);
    bestFval = bestFvalsWindow(iterationIndex);

    % Compute change in fval and individuals in last 'Window' iterations
    Window = options.StallIterLimit;
    if iteration > Window
        % The smallest fval in the window should be bestFval.
        % The largest fval in the window should be the oldest one in the
        % window. This value is at iterationIndex+1 (or 1).
        if iterationIndex == Window
            % The window runs from index 1:iterationIndex
            maxBestFvalsWindow = bestFvalsWindow(1);
        else
            % The window runs from [iterationIndex+1:end, 1:iterationIndex]
            maxBestFvalsWindow = bestFvalsWindow(iterationIndex+1);
        end
        funChange = abs(maxBestFvalsWindow-bestFval)/max(1,abs(bestFval));
    else
        funChange = Inf;
    end

    reasonToStop = '';
    exitFlag = [];

    if state.FunEval >= options.MaxFunEvals
        reasonToStop = sprintf('Optimization ended: maximum number of function evaluations exceeded.');
        exitFlag = -2;
    elseif state.Iteration >= options.MaxIter
        reasonToStop = sprintf('Optimization ended: maximum number of iterations exceeded.');
        exitFlag = -1;
    elseif bestFval <= options.ObjectiveLimit
        reasonToStop = sprintf('Optimization ended: minimum objective function limit reached.');
        exitFlag = 1;
    elseif funChange <= options.TolFun
        reasonToStop = sprintf('Optimization ended: change in the objective value less than options.TolFun.');
        exitFlag = 1;
    end

    if ~isempty(exitFlag)
        fprintf('%s\n',reasonToStop);
    end
end