clc;

clear ;

close ;

%%  Problem Definition

CostFunction = @(x) sphere(x);  % Cost Function

nVar = 5;                       % Dimension of Decision Variables

VarSize = [1,nVar];             % Matrix Size of Decision Variables

VarMin = -10;                   % Lower Bound of Decision Variables

VarMax =  10;                   % Upper Bound of Decision Variables

%%  Parameters of PSO

MaxIt = 1000;            % Maximum Number of Iterations

nPop = 50;               % Population Size 

w  = 1;                  % Inertia Coefficient

wdamp = 0.81;            % Damping Ratio of Inertia Coefficient

c1 = 2;                  % Personal Acceleration Coefficient

c2 = 2;                  % Social Acceleration Coefficient

%%  Initialization

% The Patticle Template

empty_partical.Position = [];

empty_partical.Velocity = [];

empty_partical.Cost = [];

empty_partical.Best.Position = [];

empty_partical.Best.Cost = [];

% Create Population Array

particle = repmat(empty_partical,nPop,1);

% Initialize Global Best

GlobalBest.Cost = inf;

% Iniitialize Population Members

for i=1:nPop

    % Generate Random Solution

    particle(i).Position = unifrnd(VarMin,VarMax,VarSize);

    % Initialize Velocity

    particle(i).Velocity = zeros(VarSize);

    % Evaluation

    particle(i).Cost = CostFunction(particle(i).Position);

    % Update the Personal Best

    particle(i).Best.Position = particle(i).Position;

    particle(i).Best.Cost = particle(i).Cost;

    % Update Global Best

    if particle(i).Best.Cost < GlobalBest.Cost

        GlobalBest = particle(i).Best;

    end

end

% Array to Hold Best Cost Value

BestCosts = zeros(MaxIt,1);

%%  Main Loop of PSO

for it=1:MaxIt

    for i=1:nPop

        % Update Velocity

        particle(i).Velocity =  w*particle(i).Velocity ...

          + c1*rand(VarSize).*(particle(i).Best.Position - particle(i).Position)...

          + c2*rand(VarSize).*(GlobalBest.Position - particle(i).Position);

        % Update Position

        particle(i).Position = particle(i).Position + particle(i).Velocity;

        % Evaluation

        particle(i).Cost = CostFunction( particle(i).Position);

        % Update Personal Best

        if particle(i).Cost < particle(i).Best.Cost

            particle(i).Best.Position = particle(i).Position;

            particle(i).Best.Cost = particle(i).Cost;

            % Update Global Best

            if particle(i).Best.Cost < GlobalBest.Cost

                GlobalBest = particle(i).Best;

            end

        end

    end

    % Store the Best Cost Value

    BestCosts(it) = GlobalBest.Cost;

    % Display Iteration Information

    disp(['Iteration ' num2str(it) ': Best Cost = ' num2str(BestCosts(it))]);

    % Damping Inertia Coefficient

    w = w * wdamp;

end

%%  Results

figure;

plot(BestCosts,'LineWidth',2);

semilogy(BestCosts,'LineWidth',2);

xlabel('Iterations');

ylabel('Best Cost');

grid on;

function z = sphere(x)

%%  目标函数

    z = sum(x.^2);

end