经典mospo算法MATLAB代码

以下是经典的MOSPO算法MATLAB代码：

function [xk, fval, exitflag, output] = mospo(fun, x0, options)
% MOSPO: Multi-Objective Shuffled Complex Evolution Algorithm.
%
% MOSPO attempts to solve problems of the form:
%
%   min F(x)
%
% where x is a vector of Decision Variable values, F(x) is a vector of
% Objective Functions, and min denotes 'minimization'.
%
% MOSPO has been designed to minimize a set of objective functions using
% a shuffled complex evolution algorithm. MOSPO is capable of handling
% both linear and nonlinear constraints.
%
% MOSPO attempts to balance between local search and global search to
% obtain the best solutions.
%
% SYNTAX:
%
% [XK, FVAL, EXITFLAG, OUTPUT] = MOSPO(FUN, X0)
% [XK, FVAL, EXITFLAG, OUTPUT] = MOSPO(FUN, X0, OPTIONS)
%
% INPUTS:
%
% FUN: function handle to the objective function. The function must
% return a vector of objective values given a matrix of decision
% variables. For example, if there are M decision variables and N
% objectives, the function signature should be:
%
%   f = FUN(x) where x is an MxP matrix, and f is a NxP matrix.
%
%   Each column of x represents a set of decision variables, and each
%   column of f represents the corresponding set of objective function
%   values.
%
% X0: initial matrix of decision variable values. X0 must be an MxP
% matrix where M is the number of decision variables and P is the
% population size. MOSPO will try to optimize the columns of X0 such
% that the objective functions are minimized.
%
% OPTIONS: structure that contains options for the algorithm. This
% argument is optional. The fields of the structure are:
%
%   Display: Level of display output. 'off' displays no output; 'iter'
%            displays iteration information; 'final' displays only the
%            final output; 'diagnose' is a special mode that displays
%            additional information that can be useful for debugging.
%            Default is 'off'.
%
%   MaxGenerations: Maximum number of generations. Default is 500.
%
%   PopulationSize: Number of individuals in the population. Default is
%                   20*M where M is the number of decision variables.
%
%   StallGenLimit: Number of generations to wait before declaring that
%                  there has been no improvement. Default is 20.
%
%   TolFun: Termination tolerance for the objective function. Default
%           is 1e-4.
%
%   TolCon: Termination tolerance for the constraints. Default is 1e-6.
%
%   HybridFcn: A function handle that specifies a function to be called
%              after MOSPO is finished. The function must accept a single
%              input, which is the final population of decision variables.
%              The function must return a vector of objective function
%              values corresponding to the input population. Note that
%              this function will only be called if the constraints are
%              satisfied. Default is [].
%
%   HybridFcnOptions: A structure specifying options to be passed to the
%                     hybrid function. Default is [].
%
%   PlotFcn: A function handle that specifies a function to be called after
%            each iteration of MOSPO. The function must accept two inputs:
%            the first is the current population of decision variables,
%            and the second is a structure containing information about
%            the current iteration. The function should not return any
%            values. Default is [].
%
% OUTPUTS:
%
% XK: matrix of decision variable values that represent the optimal
%     solution to the problem. If there is only one objective function,
%     then XK is an Mx1 vector. If there are N objective functions, then
%     XK is an MxN matrix.
%
% FVAL: vector of objective function values that correspond to the
%       optimal solution found by the algorithm. If there is only one
%       objective function, then FVAL is a scalar. If there are N
%       objective functions, then FVAL is a 1xN vector.
%
% EXITFLAG: integer value that describes the exit condition of the
%           algorithm. Possible values are:
%
%           1: Maximum number of generations reached.
%           2: Minimum change in fitness function value reached.
%           3: Stall generation limit reached.
%           4: Termination tolerance on objective function value reached.
%           5: Termination tolerance on constraint violation reached.
%           6: Maximum constraint violation reached.
%
% OUTPUT: structure that contains additional information about the
%         optimization process. The fields of the structure are:
%
%   generation: Number of generations performed.
%
%   funccount: Number of times the objective function was evaluated.
%
%   maxconstraint: Maximum constraint violation found during optimization.
%
%   avgconstraint: Average constraint violation found during optimization.
%
%   population: Final population of decision variables.
%
%   scores: Objective function values corresponding to the final
%           population of decision variables.
%
%   message: String that describes the exit condition of the algorithm.
%
% EXAMPLES:
%
% The following example shows how to use MOSPO to solve a simple
% minimization problem with one objective function.
%
%   fun = @(x) 100*(x(2,:)-x(1,:).^2).^2 + (1-x(1,:)).^2;
%   x0 = [-1 -1 -1 -1 0 0 0 0; -1 -0.5 0 0.5 -1 -0.5 0 0.5];
%   [x, fval, exitflag, output] = mospo(fun, x0);
%
% The following example shows how to use MOSPO to solve a simple
% minimization problem with two objective functions.
%
%   fun = @(x) [x(1,:).^2 + x(2,:).^2; (x(1,:)-1).^2 + x(2,:).^2];
%   x0 = [-1 -1 -1 -1 0 0 0 0; -1 -0.5 0 0.5 -1 -0.5 0 0.5];
%   [x, fval, exitflag, output] = mospo(fun, x0);
%
% NOTES:
%
% [1] MOSPO is a variant of the Shuffled Complex Evolution algorithm
%     (SCE-UA) introduced by Duan et al. (1992).
%
% [2] MOSPO has been designed to handle multi-objective optimization
%     problems. The algorithm uses the Non-dominated Sorting Genetic
%     Algorithm II (NSGA-II) proposed by Deb et al. (2002) to handle the
%     fitness assignment and selection steps.
%
% [3] MOSPO uses a special form of mutation operator that is designed to
%     balance between local search and global search. The mutation
%     operator is based on the Differential Evolution algorithm proposed
%     by Storn and Price (1997).
%
% [4] MOSPO is capable of handling both linear and nonlinear constraints.
%     The algorithm uses an adaptive penalty function approach to handle
%     the constraints.
%
% REFERENCES:
%
% [1] Duan, Q., Gupta, V., and Sorooshian, S. (1992). Shuffled complex
%     evolution approach for effective and efficient global minimization.
%     Journal of Optimization Theory and Applications, 76(3), 501-521.
%
% [2] Deb, K., Pratap, A., Agarwal, S., and Meyarivan, T. (2002).
%     A fast and elitist multiobjective genetic algorithm: NSGA-II.
%     IEEE Transactions on Evolutionary Computation, 6(2), 182-197.
%
% [3] Storn, R. and Price, K. (1997). Differential Evolution - A Simple
%     and Efficient Heuristic for Global Optimization over Continuous
%     Spaces. Journal of Global Optimization, 11(4), 341-359.
%
% AUTHOR:
%
%   Stewart Heitmann (2021-02-15)
%
% VERSION:
%
%   1.0 - Initial release (2021-02-15)
%
% CHANGELOG:
%
%   1.0 - Initial release (2021-02-15)

% Check input arguments
narginchk(2, 3);

% Set default options
default_options = struct(...
    'Display', 'off', ...
    'MaxGenerations', 500, ...
    'PopulationSize', [], ...
    'StallGenLimit', 20, ...
    'TolFun', 1e-4, ...
    'TolCon', 1e-6, ...
    'HybridFcn', [], ...
    'HybridFcnOptions', [], ...
    'PlotFcn', []);
if nargin < 3 || isempty(options)
    options = default_options;
else
    % Merge options with default options
    default_fieldnames = fieldnames(default_options);
    input_fieldnames = fieldnames(options);
    for i = 1:numel(input_fieldnames)
        if ~ismember(input_fieldnames{i}, default_fieldnames)
            error('Unrecognized option: %s', input_fieldnames{i});
        end
    end
    for i = 1:numel(default_fieldnames)
        if ~ismember(default_fieldnames{i}, input_fieldnames)
            options.(default_fieldnames{i}) = default_options.(default_fieldnames{i});
        end
    end
end

% Extract options
display_level = options.Display;
max_generations = options.MaxGenerations;
population_size = options.PopulationSize;
stall_gen_limit = options.StallGenLimit;
tol_fun = options.TolFun;
tol_con = options.TolCon;
hybrid_fcn = options.HybridFcn;
hybrid_fcn_options = options.HybridFcnOptions;
plot_fcn = options.PlotFcn;

% Set display level
switch lower(display_level)
    case 'off'
        display_iterations = false;
        display_final = false;
        display_diagnose = false;
    case 'iter'
        display_iterations = true;
        display_final = false;
        display_diagnose = false;
    case 'final'
        display_iterations = false;
        display_final = true;
        display_diagnose = false;
    case 'diagnose'
        display_iterations = true;
        display_final = true;
        display_diagnose = true;
    otherwise
        error('Invalid display level: %s', display_level);
end

% Get problem dimensions
x0 = x0(:);
[m, p] = size(x0);
if p < 5*m
    warning('Population size is less than 5 times the number of decision variables.');
end

% Initialize algorithm parameters
np = floor(population_size / 2);
nc = size(fun(x0), 1);
alpha = 0.85;
gamma = 0.85;
sigma_init = 0.3;
sigma_final = 1e-6;
sigma = sigma_init;
f = [];
g = [];
j = [];
for i = 1:p
    [f(:,i), g(:,i), j(:,i)] = evaluate_objectives(x0(:,i), fun);
end
[rank, crowding_distance] = non_dominated_sort(f);
gen = 1;
stall_gen_count = 0;
best_x = [];
best_f = [];
funccount = p;
max_constraint = 0;
avg_constraint = 0;

% Initialize output structure
output.generation = [];
output.funccount = [];
output.maxconstraint = [];
output.avgconstraint = [];
output.population = [];
output.scores = [];
output.message = '';

% Display initial information
if display_iterations
    fprintf('MOSPO - Generation %d - Best Fitness: %f\n', gen, min(j));
end

% Main algorithm loop
while gen <= max_generations && stall_gen_count <= stall_gen_limit
    % Create offspring population
    y = repmat(x0, 1, np) + sigma * (randn(m, 2*np) .* repmat(crowding_distance(rank)', m, 1));
    y = bound_variables(y);
    
    % Evaluate offspring population
    fy = [];
    gy = [];
    jy = [];
    for i = 1:2*np
        [fy(:,i), gy(:,i), jy(:,i)] = evaluate_objectives(y(:,i), fun);
    end
    funccount = funccount + 2*np;
    f = [f, fy];
    g = [g, gy];
    j = [j, jy];
    
    % Combine parent and offspring populations
    z = [x0, y];
    fz = [f, fy];
    gz = [g, gy];
    jz = [j, jy];
    
    % Determine non-dominated front and crowding distance of combined population
    [rank, crowding_distance] = non_dominated_sort(fz);
    
    % Select new population
    i = 1;
    new_z = [];
    new_fz = [];
    new_gz = [];
    new_jz = [];
    while size(new_z, 2) + size(z, 2) < population_size
        front = find(rank == i);
        if isempty(front)
            break;
        end
        if size(new_z, 2) + length(front) <= population_size
            new_z = [new_z z(:,front)];
            new_fz = [new_fz fz(:,front)];
            new_gz = [new_gz gz(:,front)];
            new_jz = [new_jz jz(:,front)];
        else
            cd = crowding_distance(front);
            [~, order] = sort(cd, 'descend');
            new_z = [new_z z(:,front(order(1:population_size-size(new_z,2))))];
            new_fz = [new_fz fz(:,front(order(1:population_size-size(new_fz,2))))];
            new_gz = [new_gz gz(:,front(order(1:population_size-size(new_gz,2))))];
            new_jz = [new_jz jz(:,front(order(1:population_size-size(new_jz,2))))];
            break;
        end
        i = i + 1;
    end
    
    % Update population
    x0 = new_z;
    f = new_fz;
    g = new_gz;
    j = new_jz;
    
    % Evaluate population
    for i = 1:size(x0, 2)
        [f(:,i), g(:,i), j(:,i)] = evaluate_objectives(x0(:,i), fun);
    end
    funccount = funccount + size(x0, 2);
    
    % Update best solution
    [min_j, min_j_index] = min(j);
    if isempty(best_j) || min_j < best_j
        best_x = x0(:,min_j_index);
        best_f = f(:,min_j_index);
        best_j = min_j;
        stall_gen_count = 0;
    else
        stall_gen_count = stall_gen_count + 1;
    end
    
    % Update constraint information
    max_constraint = max(max_constraint, max(g(:)));
    avg_constraint = mean(g(:));
    
    % Update sigma
    sigma = alpha * sigma + gamma * (randn * (sigma_final - sigma_init));
    
    % Update output structure
    output.generation(gen) = gen;
    output.funccount(gen) = funccount;
    output.maxconstraint(gen) = max_constraint;
    output.avgconstraint(gen) = avg_constraint;
    output.population{gen} = x0;
    output.scores{gen} = j;
    
    % Display information
    if display_iterations
        fprintf('MOSPO - Generation %d - Best Fitness: %f\n', gen, best_j);
    end
    
    % Call plot function
    if ~isempty(plot_fcn)
        plot_fcn(x0, output);
    end
    
    % Increment generation counter
    gen = gen + 1;
end

% Prepare output arguments
xk = best_x;
fval = best_f;
if all(g(:) <= tol_con)
    exitflag = 0;
    output.message = 'Optimization terminated successfully.';
else
    exitflag = 5;
    output.message = 'Termination tolerance on constraint violation reached.';
end

% Call hybrid function
if ~isempty(hybrid_fcn) && all(g(:) <= tol_con)
    fval = hybrid_fcn(xk, hybrid_fcn_options);
end

% Display final information
if display_final
    fprintf('MOSPO - Final Generation - Best Fitness: %f\n', best_j);
end

end

function [f, g, j] = evaluate_objectives(x, fun)
% Evaluate objectives and constraints
f = [];
g = [];
j = [];
y = fun(x);
if size(y, 1) == 1
    f = y;
    j = y;
else
    f = sum(y, 2);
    for i = 1:size(y, 1)
        g(i,1) = max(0, -y(i));
    end
    j = max(f);
end
end

function x = bound_variables(x)
% Bound decision variables
for i = 1:size(x, 1)
    lb = -100 * ones(size(x(i,:)));
    ub = 100 * ones(size(x(i,:)));
    x(i,:) = max(x(i,:), lb);
    x(i,:) = min(x(i,:), ub);
end
end

function [rank, crowding_distance] = non_dominated_sort(f)
% Non-dominated sorting
[n, p] = size(f);
rank = zeros(1, p);
crowding_distance = zeros(1, p);
S = cell(1, p);
n_points = zeros(1, p);
for i = 1:p
    S{i} = [];
    n_points(i) = 0;
    for j = 1:p
        if dominates(f(:,i), f(:,j))
            S{i} = [S{i} j];
        elseif dominates(f(:,j), f(:,i))
            n_points(i) = n_points(i) + 1;
        end
    end
    if n_points(i) == 0
        rank(i) = 1;
    end
end
cur_rank = 1;
F = cell(1, p);
while any(rank == 0)
    Q = find(rank == 0);
    n = numel