以下是经典的MOSPO算法MATLAB代码：
```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
扫一扫
1260
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
