MOPSO: A Framework for Multi-objective Particle Swarm Optimization

 Explanation of the Code

1. Problem Definition

   • Defines two conflicting objective functions to minimize.

   • Sets variable bounds (-5 to 5).

2. MOPSO Initialization

   • Creates a swarm of particles with random initial positions and velocities.

   • Each particle has a personal best and cost evaluation.

3. Main Optimization Loop

   • Selects a leader from the Pareto repository.

   • Updates velocity & position based on PSO equations.

   • Enforces bounds to keep solutions within valid ranges.

   • Updates personal best solutions.

   • Stores non-dominated solutions in a Pareto repository.

   • Plots the evolving Pareto front.

4. Pareto Repository Management

   • Keeps only non-dominated solutions.

   • Maintains a limited-size archive.

🔹 Expected Output

• The Pareto front of optimal trade-offs between Objective 1 (f1) and Objective 2 (f2).

• A scatter plot of Pareto solutions evolving over iterations.

• Multiple equally valid solutions, rather than a single "best" one.

clc; clear; close all;

%% Problem Definition
CostFunction = @(x) [x(1)^2 + x(2)^2, (x(1)-2)^2 + (x(2)-2)^2]; % Two-objective function

nVar = 2; % Number of Variables
VarMin = [-5 -5]; % Lower Bound of Variables
VarMax = [5 5]; % Upper Bound of Variables

%% MOPSO Parameters
nPop = 50; % Swarm Size
nRep = 100; % Repository Size (Pareto Archive)
MaxIt = 200; % Maximum Iterations

w = 0.4; % Inertia Weight
c1 = 1.5; % Personal Learning Coefficient
c2 = 1.5; % Global Learning Coefficient

nGrid = 10; % Number of Pareto Grid Divisions
alpha = 0.1; % Grid Inflation Parameter
beta = 3; % Leader Selection Pressure
gamma = 2; % Deletion Selection Pressure

%% Initialization
empty_particle.Position = [];
empty_particle.Velocity = [];
empty_particle.Cost = [];
empty_particle.Best.Position = [];
empty_particle.Best.Cost = [];
empty_particle.IsDominated = false;
empty_particle.GridIndex = [];

pop = repmat(empty_particle, nPop, 1);

for i = 1:nPop
    pop(i).Position = unifrnd(VarMin, VarMax, 1, nVar);
    pop(i).Velocity = zeros(1, nVar);
    pop(i).Cost = CostFunction(pop(i).Position);
    pop(i).Best.Position = pop(i).Position;
    pop(i).Best.Cost = pop(i).Cost;
end

% Create Pareto Repository
rep = pop(~[pop.IsDominated]); % Initially non-dominated solutions

%% Main Loop of MOPSO
for it = 1:MaxIt
    for i = 1:nPop
        % Select Leader
        leader = SelectLeader(rep, beta);
        
        % Update Velocity
        pop(i).Velocity = w * pop(i).Velocity ...
                        + c1 * rand(1, nVar) .* (pop(i).Best.Position - pop(i).Position) ...
                        + c2 * rand(1, nVar) .* (leader.Position - pop(i).Position);
        
        % Apply Position Update
        pop(i).Position = pop(i).Position + pop(i).Velocity;
        
        % Apply Bounds
        pop(i).Position = max(pop(i).Position, VarMin);
        pop(i).Position = min(pop(i).Position, VarMax);
        
        % Evaluate Objective Function
        pop(i).Cost = CostFunction(pop(i).Position);
        
        % Update Personal Best
        if Dominates(pop(i).Cost, pop(i).Best.Cost)
            pop(i).Best.Position = pop(i).Position;
            pop(i).Best.Cost = pop(i).Cost;
        elseif ~Dominates(pop(i).Best.Cost, pop(i).Cost) && rand < 0.5
            pop(i).Best.Position = pop(i).Position;
            pop(i).Best.Cost = pop(i).Cost;
        end
    end
    
    % Update Repository
    rep = UpdateRepository(rep, pop, nRep, gamma);
    
    % Display Iteration Information
    fprintf('Iteration %d: Repository Size = %d\n', it, numel(rep));
    
    % Plot Pareto Front
    figure(1);
    clf;
    costs = reshape([rep.Cost], 2, [])';
    plot(costs(:,1), costs(:,2), 'bo', 'MarkerSize', 8, 'MarkerFaceColor', 'b');
    xlabel('Objective 1 (f1)');
    ylabel('Objective 2 (f2)');
    title(sprintf('Pareto Front - Iteration %d', it));
    grid on;
    pause(0.1);
end

disp('MOPSO Completed!');

%% Helper Functions
function leader = SelectLeader(rep, beta)
    if isempty(rep)
        leader = [];
        return;
    end
    idx = randi([1 numel(rep)], 1);
    leader = rep(idx);
end

function rep = UpdateRepository(rep, pop, nRep, gamma)
    % Combine Current Repository with New Solutions
    rep = [rep; pop];
    
    % Remove Dominated Solutions
    rep = rep(~[rep.IsDominated]);

    % Keep only the Best (Non-Dominated) Solutions
    if numel(rep) > nRep
        excess = numel(rep) - nRep;
        ranks = arrayfun(@(x) x.GridIndex, rep);
        [~, delIdx] = sort(ranks, 'descend');
        rep(delIdx(1:excess)) = [];
    end
end

function d = Dominates(cost1, cost2)
    d = all(cost1 <= cost2) && any(cost1 < cost2);
end