Can't extract variables that give the best cost function from PSO code
5 views (last 30 days)
Show older comments
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
3 Comments
Sam Chak
on 25 Aug 2023
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.
Answers (1)
Sam Chak
on 25 Aug 2023
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.
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
0 Comments
See Also
Categories
Find more on Particle Swarm in Help Center and File Exchange
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!