一、运行结果
二、多目标粒子群算法（MOPSO）函数寻优Matlab代码
clc;
clear;
close all;
%% Problem Definition
CostFunction = @(x) ZDT(x); % Cost Function
nVar = 5; % Number of Decision Variables
VarSize = [1 nVar]; % Size of Decision Variables Matrix
VarMin = 0; % Lower Bound of Variables
VarMax = 1; % Upper Bound of Variables
%% MOPSO Parameters
MaxIt = 200; % Maximum Number of Iterations
nPop = 200; % Population Size
nRep = 100; % Repository Size
w = 0.5; % Inertia Weight
wdamp = 0.99; % Intertia Weight Damping Rate
c1 = 1; % Personal Learning Coefficient
c2 = 2; % Global Learning Coefficient
nGrid = 7; % Number of Grids per Dimension
alpha = 0.1; % Inflation Rate
beta = 2; % Leader Selection Pressure
gamma = 2; % Deletion Selection Pressure
mu = 0.1; % Mutation Rate
%% Initialization
empty_particle.Position = [];
empty_particle.Velocity = [];
empty_particle.Cost = [];
empty_particle.Best.Position = [];
empty_particle.Best.Cost = [];
empty_particle.IsDominated = [];
empty_particle.GridIndex = [];
empty_particle.GridSubIndex = [];
pop = repmat(empty_particle, nPop, 1);
for i = 1:nPop
pop(i).Position = unifrnd(VarMin, VarMax, VarSize);
pop(i).Velocity = zeros(VarSize);
pop(i).Cost = CostFunction(pop(i).Position);
% Update Personal Best
pop(i).Best.Position = pop(i).Position;
pop(i).Best.Cost = pop(i).Cost;
end
% Determine Domination
pop = DetermineDomination(pop);
rep = pop(~[pop.IsDominated]);
Grid = CreateGrid(rep, nGrid, alpha);
for i = 1:numel(rep)
rep(i) = FindGridIndex(rep(i), Grid);
end
%% MOPSO Main Loop
for it = 1:MaxIt
for i = 1:nPop
leader = SelectLeader(rep, beta);
pop(i).Velocity = w*pop(i).Velocity ...
+c1*rand(VarSize).*(pop(i).Best.Position-pop(i).Position) ...
+c2*rand(VarSize).*(leader.Position-pop(i).Position);
pop(i).Position = pop(i).Position + pop(i).Velocity;
pop(i).Position = max(pop(i).Position, VarMin);
pop(i).Position = min(pop(i).Position, VarMax);
pop(i).Cost = CostFunction(pop(i).Position);
if Dominates(pop(i), pop(i).Best)
pop(i).Best.Position = pop(i).Position;
pop(i).Best.Cost = pop(i).Cost;
elseif Dominates(pop(i).Best, pop(i))
% Do Nothing
else
if rand<0.5
pop(i).Best.Position = pop(i).Position;
pop(i).Best.Cost = pop(i).Cost;
end
end
end
% Add Non-Dominated Particles to REPOSITORY
rep = [rep
pop(~[pop.IsDominated])]; %#ok
% Determine Domination of New Resository Members
rep = DetermineDomination(rep);
% Keep only Non-Dminated Memebrs in the Repository
rep = rep(~[rep.IsDominated]);
% Update Grid
Grid = CreateGrid(rep, nGrid, alpha);
% Update Grid Indices
for i = 1:numel(rep)
rep(i) = FindGridIndex(rep(i), Grid);
end
% Check if Repository is Full
if numel(rep)>nRep
Extra = numel(rep)-nRep;
for e = 1:Extra
rep = DeleteOneRepMemebr(rep, gamma);
end
end
% Plot Costs
figure(1);
PlotCosts(pop, rep);
pause(0.01);
% Show Iteration Information
disp(['Iteration ' num2str(it) ': Number of Rep Members = ' num2str(numel(rep))]);
% Damping Inertia Weight
w = w*wdamp;
end
function particle = FindGridIndex(particle, Grid)
nObj = numel(particle.Cost);
nGrid = numel(Grid(1).LB);
particle.GridSubIndex = zeros(1, nObj);
for j = 1:nObj
particle.GridSubIndex(j) = ...
find(particle.Cost(j)<Grid(j).UB, 1, 'first');
end
particle.GridIndex = particle.GridSubIndex(1);
for j = 2:nObj
particle.GridIndex = particle.GridIndex-1;
particle.GridIndex = nGrid*particle.GridIndex;
particle.GridIndex = particle.GridIndex+particle.GridSubIndex(j);
end
end
function b = Dominates(x, y)
if isstruct(x)
x = x.Cost;
end
if isstruct(y)
y = y.Cost;
end
b = all(x <= y) && any(x<y);
end
function pop = DetermineDomination(pop)
nPop = numel(pop);
for i = 1:nPop
pop(i).IsDominated = false;
end
for i = 1:nPop-1
for j = i+1:nPop
if Dominates(pop(i), pop(j))
pop(j).IsDominated = true;
end
if Dominates(pop(j), pop(i))
pop(i).IsDominated = true;
end
end
end
end
function rep = DeleteOneRepMemebr(rep, gamma)
% Grid Index of All Repository Members
GI = [rep.GridIndex];
% Occupied Cells
OC = unique(GI);
% Number of Particles in Occupied Cells
N = zeros(size(OC));
for k = 1:numel(OC)
N(k) = numel(find(GI == OC(k)));
end
% Selection Probabilities
P = exp(gamma*N);
P = P/sum(P);
% Selected Cell Index
sci = RouletteWheelSelection(P);
% Selected Cell
sc = OC(sci);
% Selected Cell Members
SCM = find(GI == sc);
% Selected Member Index
smi = randi([1 numel(SCM)]);
% Selected Member
sm = SCM(smi);
% Delete Selected Member
rep(sm) = [];
end
function Grid = CreateGrid(pop, nGrid, alpha)
c = [pop.Cost];
cmin = min(c, [], 2);
cmax = max(c, [], 2);
dc = cmax-cmin;
cmin = cmin-alpha*dc;
cmax = cmax+alpha*dc;
nObj = size(c, 1);
empty_grid.LB = [];
empty_grid.UB = [];
Grid = repmat(empty_grid, nObj, 1);
for j = 1:nObj
cj = linspace(cmin(j), cmax(j), nGrid+1);
Grid(j).LB = [-inf cj];
Grid(j).UB = [cj +inf];
end
end
function leader = SelectLeader(rep, beta)
% Grid Index of All Repository Members
GI = [rep.GridIndex];
% Occupied Cells
OC = unique(GI);
% Number of Particles in Occupied Cells
N = zeros(size(OC));
for k = 1:numel(OC)
N(k) = numel(find(GI == OC(k)));
end
% Selection Probabilities
P = exp(-beta*N);
P = P/sum(P);
% Selected Cell Index
sci = RouletteWheelSelection(P);
% Selected Cell
sc = OC(sci);
% Selected Cell Members
SCM = find(GI == sc);
% Selected Member Index
smi = randi([1 numel(SCM)]);
% Selected Member
sm = SCM(smi);
% Leader
leader = rep(sm);
end
function i = RouletteWheelSelection(P)
r = rand;
C = cumsum(P);
i = find(r <= C, 1, 'first');
end
function z = ZDT(x)
n = numel(x);
f1 = x(1);
g = 1+9/(n-1)*sum(x(2:end));
h = 1-sqrt(f1/g);
f2 = g*h;
z = [f1
f2];
end
function PlotCosts(pop, rep)
pop_costs = [pop.Cost];
plot(pop_costs(1, :), pop_costs(2, :), 'ko');
hold on;
rep_costs = [rep.Cost];
plot(rep_costs(1, :), rep_costs(2, :), 'r*');
xlabel('1^{st} Objective');
ylabel('2^{nd} Objective');
grid on;
hold off;
end