Can't extract variables that give the best cost function from PSO code

Good day. I am currently using a particle swarm algorithm to optimize PID gains. The PSO function is attached. I want to get the PID gains that result in the best cost function value. This value is stored in BestVar. I am getting something out but it does not give me the best value displayed in BestCosts. I'm sure BestCosts displays the correct variables as I've also used the built-in particle swarm function. I can't figure out where I am going wrong. Please assist.

function output=PSO(problem, params)

    %% Problem definition
    
    costfunction = problem.costfunction;   % Cost Function
    
    nVar = problem.nVar;             % Number of unknown (decision) variables;
    varSize = [1 nVar];              % Matrix size of decision variables
    varMin = problem.varMin;         % lower bound of decision variables
    varMax = problem.varMax;         % upper bound of decision variables
    
    %% Parameters of PSO
    
    maxIt = params.maxIt;                    % maximum number of iterations
    nPop = params.nPop;                       % population size/ swarm size
    w = params.w;                           % inertial coefficient (standard value)
    wdamp = params.wdamp;                     % Damping ratio of inertia weight. Reducing inertia coefficient makes performanc of algorithm better
    c1 = params.c1;                          % personal acceleration (learning) coefficient
    c2 = params.c2;                          % social/global acceleration coefficient

    % the flag for showing iteration information
    %showIterInfo = params.showIterInfo; 
    
    maxVelocity = 0.2*(varMax - varMin);
    minVelocity = -maxVelocity;

    %% Initialization
    
    % The particle template
    empty_particle.Position = [];    %  Position of an emplty (initialization ??) particle
    empty_particle.Velocity = [];    % Velocity of empty particle
    empty_particle.Cost = [];        % Particle cost value
    empty_particle.Best.Position = [];  % Position of the best particle
    empty_particle.Best.Cost = [];      % Cost of the best particle
    
    %repeat above for the whole population, i.e Create population array
    
    particle = repmat(empty_particle, nPop, 1); % contains fields Position, Velocity, Cost, Best Position Best Cost
    
    % Initialize global best
    GlobalBest.Cost = inf;  % for minimization, global best is +inf (worst value). For maximization, global best is -inf 
    
    % Initialize population members
    for i = 1:nPop
        
        % Generate random solution
        particle(i).Position = unifrnd(varMin,varMax,varSize);  % uniform distribution
    
        % 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 on each iteration
    BestCosts = zeros(maxIt,1);
    BestVar = zeros(maxIt,nVar);
    
    %% 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);

            % Apply velocity limits
            particle(i).Velocity = max(particle(i).Velocity, minVelocity);
            particle(i).Velocity = min(particle(i).Velocity, maxVelocity);

            % Update position
            particle(i).Position = particle(i).Position + particle(i).Velocity;

            % Apply lower and upper bound limits
            particle(i).Position = max(particle(i).Position, varMin);
            particle(i).Position = min(particle(i).Position, varMin);
            
            % 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;
        BestVar(it,:) = GlobalBest.Position;
    
        % Display iteration information
        if params.dispIterInfo==true
            disp(['Iteration' num2str(it) ': Best Cost =' num2str(BestCosts(it))])
            hold on
            plot(it,BestCosts(it), 'b.');
            xlim([0 maxIt]);
            title(['Best Cost =',num2str(BestCosts(it))]);
            xlabel('Iteration');
            ylabel('Best Cost');
        end
    
        % Damping inertia coefficient
        w = w*wdamp;
    
    end

    output.pop = particle;
    output.GlobalBest = GlobalBest;
    output.BestVar=BestVar;
    output.BestCosts= BestCosts;

end

[Sam Chak]

Can you show the cost function so that we can duplicate the same PID gains like you do?

[Phumelele Magagula]

Hi Sam.

The function draws from a simulink model, so I'm afraid it won't be of much use. I'll post it anyway

function cost=optimization_PID(k)
simulationdata=sim('LOS_Controller_09_08');
b=simulationdata;
t=b.tout;
J_err_out=b.J_err;
J_err=max(J_err_out(:,end,end));

assignin('base','k',k);
sim("LOS_Controller_09_08");
cost=J_err(length(J_err));
end

[Sam Chak]

Hi @Phumelele Magagula

I have attached a Simulink file as a toy problem for you to test your code. Please show exactly how you run the PSO optimization code. Your code requires:

problem.costfunction

Copy and paste the results.

Answer 1

[Sam Chak]

Hi @Phumelele Magagula

I've made some modifications to the PSO code in order to test the toy problem, as certain crucial pieces of information (problem, params) were not provided. I can assure you that the code is capable of identifying specific suboptimal values for the PID gains and the cost. However, please note that you'll need to provide your cost function in order to test your own system effectively.

problem = @PIDcostfun;

% output = PSO(problem);
% fNames = fieldnames(output);
% for n  = 1:length(fNames)
%     disp(output.(fNames{n}))
% end

GlobalBest = PSO(problem)   % True optimal gain values are [1 1 1] and Best Cost is 3.

GlobalBest = struct with fields:
    Position: [1.0744 1.0152 1.0590]
        Cost: 3.0018

function GlobalBest = PSO(problem)
% function output = PSO(problem)

    %% Problem definition
    
    costfunction = problem;   % Cost Function
    
    nVar    = 3;                % Number of unknown (decision) variables;
    varSize = [1 nVar];         % Matrix size of decision variables
    varMin  = [0.5 0.5 0.5];    % lower bound of decision variables
    varMax  = [1.5 1.5 1.5];    % upper bound of decision variables
    
    %% Parameters of PSO
    
    maxIt = 200;        % maximum number of iterations
    nPop  = 100;        % population size/ swarm size
    w     = 1;          % inertial coefficient (standard value)
    wdamp = 0.99;       % Damping ratio of inertia weight. Reducing inertia coefficient makes performanc of algorithm better
    c1    = 1.50;       % personal acceleration (learning) coefficient
    c2    = 2.00;       % social/global acceleration coefficient

    % the flag for showing iteration information
    %showIterInfo = params.showIterInfo; 
    
    maxVelocity = 0.2*(varMax - varMin);
    minVelocity = -maxVelocity;

    %% Initialization
    
    % The particle template
    empty_particle.Position = [];    %  Position of an emplty (initialization ??) particle
    empty_particle.Velocity = [];    % Velocity of empty particle
    empty_particle.Cost = [];        % Particle cost value
    empty_particle.Best.Position = [];  % Position of the best particle
    empty_particle.Best.Cost = [];      % Cost of the best particle
    
    %repeat above for the whole population, i.e Create population array
    
    particle = repmat(empty_particle, nPop, 1); % contains fields Position, Velocity, Cost, Best Position Best Cost
    
    % Initialize global best
    GlobalBest.Cost = inf;  % for minimization, global best is +inf (worst value). For maximization, global best is -inf 
    
    % Initialize population members
    for i = 1:nPop
        
        % Generate random solution
        particle(i).Position = unifrnd(varMin,varMax,varSize);  % uniform distribution
    
        % 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 on each iteration
    BestCosts = zeros(maxIt,1);
    BestVar = zeros(maxIt,nVar);
    
    %% 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);

            % Apply velocity limits
            particle(i).Velocity = max(particle(i).Velocity, minVelocity);
            particle(i).Velocity = min(particle(i).Velocity, maxVelocity);

            % Update position
            particle(i).Position = particle(i).Position + particle(i).Velocity;

            % Apply lower and upper bound limits
            particle(i).Position = max(particle(i).Position, varMin);
            particle(i).Position = min(particle(i).Position, varMin);
            
            % 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;
        BestVar(it,:) = GlobalBest.Position;
    
        % Display iteration information
        % if params.dispIterInfo==true
        %     disp(['Iteration' num2str(it) ': Best Cost =' num2str(BestCosts(it))])
        %     hold on
        %     plot(it,BestCosts(it), 'b.');
        %     xlim([0 maxIt]);
        %     title(['Best Cost =',num2str(BestCosts(it))]);
        %     xlabel('Iteration');
        %     ylabel('Best Cost');
        % end
    
        % Damping inertia coefficient
        w = w*wdamp;
    
    end

    % output.pop = particle;
    % output.GlobalBest = GlobalBest;
    % output.BestVar=BestVar;
    % output.BestCosts= BestCosts;

end