多目标优化算法之多目标黏菌算法MOSMA（附Matlab免费代码）-优快云博客

本文链接：https://blog.youkuaiyun.com/weixin_44028734/article/details/145110542

引言

Multi-Objective Slime Mould Algorithm (MOSMA) 是一种新兴的多目标优化算法，受到自然界中黏菌行为的启发。黏菌在寻找食物时表现出的动态行为和优化路径选择能力，为设计高效的优化算法提供了灵感。MOSMA 通过模拟黏菌的生命周期中的营养、活跃和动态阶段，实现了对多目标优化问题的求解。

1. 算法原理

MOSMA 的核心思想是通过模拟黏菌在寻找食物时的行为模式，包括接近食物、包围食物和消化食物等过程，来实现优化。具体来说，算法通过以下步骤进行：

初始化：生成一组随机的候选解，每个解代表一个黏菌的位置。
适应度评估：计算每个候选解的适应度值，通常使用目标函数的值来评估。
排序：根据适应度值对候选解进行排序，通常使用非支配排序方法。
权重更新：根据排序结果更新每个候选解的权重，权重较高的解表示更优的解。
位置更新：根据权重和随机选择的其他候选解，更新每个候选解的位置。

2. 应用场景

机械工程设计：用于解决复杂的机械工程优化问题，如电机设计和自动电压调节系统
图像处理：用于图像分割和特征选择，提高图像处理的效率和准确性
电力系统：用于电力系统的谐波滤波器设计，优化电力系统的性能

从代距（GD）、倒代距（IGD）、最大传播（MS）、间隔和运行时间等性能指标对MOSMA进行了性能评估。仿真结果表明，该算法在求解线性、非线性、连续和离散Pareto最优前沿等多目标问题方面具有优势。

参考文献

Matlab代码下载

微信搜索并关注-优化算法侠（英文名：Swarm-Opti），或扫描下方二维码关注，以算法名字搜索历史文章即可下载。


% 请关注微信公众号：优化算法侠，发现更多精彩

function f = MOSMA(dim,M,lb,ub,N,Max_iter)          
X = zeros(N,dim);
Sol = zeros(N,dim);
weight = ones(N,dim);%fitness weight of each slime mold
%% Initialize the population
for i=1:N
   x(i,:)=lb+(ub-lb).*rand(1,dim); 
   f(i,1:M) = evaluate_objective(x(i,:), M);
end
new_Sol=[x f]; 
new_Sol = solutions_sorting(new_Sol, M, dim);
for i = 1 : Max_iter 
[SmellOrder,SmellIndex] = sort(Sol);  
worstFitness = SmellOrder(N);
bestFitness = SmellOrder(1);
S=bestFitness-worstFitness+eps;  % plus eps to avoid denominator zero
        for k=1:N
            if k<=(N/2)  
                weight(SmellIndex(k),:) = 1+rand()*log10((bestFitness-SmellOrder(k))/(S)+1);
            else
                weight(SmellIndex(k),:) = 1-rand()*log10((bestFitness-SmellOrder(k))/(S)+1);
            end
        end     
a = atanh(-(i/Max_iter)+1);   
b = 1-i/Max_iter;
    for j=1:N 
        best=(new_Sol(j,1:dim) - new_Sol(1,(1:dim)));
        if rand<0.03    
            X(j,:) = (ub-lb).*rand+lb;
        else
            p =tanh(abs(f(j)-best));  
            vb = unifrnd(-a,a,1,dim); 
            vc = unifrnd(-b,b,1,dim);        
                r = rand();
                A = randi([1,N]);  
                B = randi([1,N]);
                if r<p    
                    X(j,:) = best+ vb.*(weight(j,:).*X(A,:)-X(B,:));
                else
                    X(j,:) = best+ vc.*(weight(j,:).*X(A,:)-X(B,:));
                end  
        end
        Sol(j,1:dim) = X(j,1:dim);       
        Flag4ub=Sol(j,1:dim)>ub;
        Flag4lb=Sol(j,1:dim)<lb;
        Sol(j,1:dim)=(Sol(j,1:dim).*(~(Flag4ub+Flag4lb)))+ub.*Flag4ub+lb.*Flag4lb;  
        %% Evalute the fitness/function values of the new population
        Sol(j, dim+1:M+dim) = evaluate_objective(Sol(j,1:dim),M);
        if Sol(j,dim+1:dim+M) <= new_Sol(1,(dim+1:dim+M)) 
           new_Sol(1,1:(dim+M)) = Sol(j,1:(dim+M));  
        end
    end    
%% ! Very important to combine old and new bats !
   Sort_bats(1:N,:) = new_Sol;
   Sort_bats((N + 1):(2*N), 1:M+dim) = Sol;
%% Non-dominated sorting process (a separate function/subroutine)
   Sorted_bats = solutions_sorting(Sort_bats, M, dim); 
%% Select npop solutions among a combined population of 2*npop solutions  
    new_Sol = cleanup_batspop(Sorted_bats, M, dim, N);  
end
f=new_Sol;

end



%% Clean up the populations (both old and new) to give a new population
% This cleanup here is similar to the Non-dominated Sorting Genetic
% Algorithm (NSGA-II) by K. Deb et al. (2002), which can be applied to 
% any cleanup of 2*npop solutions to form a set of npop solutions.
function new_bats = cleanup_batspop(bats, m, ndim, npop)
% The input population to this part has twice (ntwice) of the needed 
% population size (npop). Thus, selection is done based on ranking and 
% crowding distances, calculated from the non-dominated sorting
ntwice= size(bats,1);
% Ranking is stored in column Krank
Krank=m+ndim+1;
% Sort the population of size 2*npop according to their ranks
[~,Index] = sort(bats(:,Krank)); sorted_bats=bats(Index,:);
% Get the maximum rank among the population
RankMax=max(bats(:,Krank)); 

%% Main loop for selecting solutions based on ranks and crowding distances
K = 0;  % Initialization for the rank counter 
% Loop over all ranks in the population
for i =1:RankMax,  
    % Obtain the current rank i from sorted solutions
    RankSol = max(find(sorted_bats(:, Krank) == i));
    % In the new bats/solutions, there can be npop solutions to fill
    if RankSol<npop,
       new_bats(K+1:RankSol,:)=sorted_bats(K+1:RankSol,:);
    end 
    % If the population after addition is large than npop, re-arrangement
    % or selection is carried out
    if RankSol>=npop
        % Sort/Select the solutions with the current rank 
        candidate_bats = sorted_bats(K + 1 : RankSol, :);
        [~,tmp_Rank]=sort(candidate_bats(:,Krank+1),'descend');
        % Fill the rest (npop-K) bats/solutions up to npop solutions 
        for j = 1:(npop-K), 
            new_bats(K+j,:)=candidate_bats(tmp_Rank(j),:);
        end
    end
    % Record and update the current rank after adding new bats 
    K = RankSol;
end
end



%% The sorting of nondomninated solutions from a population of 2*npop     %
%% (New solutions+old population) to form a population of npop solutions  %
function sorted_x = solutions_sorting(x, m, ndim)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Inputs and outputs are the extended solutions x with a dimension of    %
%% npop by (ndim+m+2). The objective values are already included in x.    %
% More specifically, the first ndim columns are the actual solutions or 
% variable values (1:ndim), followed by the m columns of objective values. 
% Then, the next column (i.e.,ndim+m+1) corresponds to the ranks, whereas 
% the final column (i.e., ndim+m+2) records the crowd distances. 
% ----------------------------------------------------------------------- %
%% Get the parameters from the inputs such as the size of input population
npop=size(x,1);    % Population size
frontRank=1;       % Pareto frontRank (counter) initialization
Rcol=ndim+m+1;     % Store the ranks in the column Rcol=ndim+m+1
% Define the Parato Front as a class (PF) and initilization of xSol
PF(frontRank).R=[];   xSol=[];
%% The main non-dominated sorting starts here             %%%%%%%%%%%%%%%%% 
for i = 1:npop, 
    % Set the number (initially, 0) of solutions dominating this solution
    xSol(i).n=0;
    % Find all the solutions (that dominated by this solution)
    xSol(i).q=[];
    % Sorting into 3 categories: better (minimization), equal & otherwise
    for j=1:npop,
        % Definte 3 counters for 3 categories
        ns_categ_1=0; ns_categ_2=0; ns_categ_3=0;
        for k=1:m,  % for all m objectives
            % Update the counters for 3 different categories
            if (x(i,ndim+k) < x(j,ndim+k)),      % better/non-dominated
                ns_categ_1=ns_categ_1+1;
            elseif (x(i,ndim+k)==x(j,ndim+k)),   % equal
                ns_categ_2=ns_categ_2+1;
            else                                 % dominated
                ns_categ_3=ns_categ_3+1;
            end
        end % end of k
        % Update the solutions in their class
        if ns_categ_1==0 && ns_categ_2 ~= m
            xSol(i).n=xSol(i).n+1;
        elseif ns_categ_3 == 0 && ns_categ_2 ~= m
            xSol(i).q=[xSol(i).q j];
        end
    end % end of j   
    %% Record/Udpate the Pareto Front
    if xSol(i).n==0, 
        x(i,Rcol)=1;   % Update the Rank #1 (i.e., the Pareto Front)
        PF(frontRank).R = [PF(frontRank).R i];
    end
end % end of i=1:npop (The first round full rank-sorting process)

% Update the rest frontRanks (close, but not on the Pareto Front)
while ~isempty(PF(frontRank).R),
    nonPF=[];    % Intialization the set
    N=length(PF(frontRank).R);
for i=1 :N, 
   % Get the solution/list 
   Sol_tmp_q=xSol(PF(frontRank).R(i)).q; 
   % If not empty, update 
   if ~isempty(xSol(Sol_tmp_q))
       for j = 1:length(Sol_tmp_q),
         % Get the solutions dominated by the current solution    
          Sol_tmp_qj=xSol(PF(frontRank).R(i)).q(j);   
          xSol(Sol_tmp_qj).n=xSol(Sol_tmp_qj).n-1;
          if xSol(Sol_tmp_qj).n==0
             x(Sol_tmp_qj, Rcol)=frontRank + 1;
             nonPF = [nonPF Sol_tmp_qj];
          end
       end % end of j
   end
end  % end of i
   frontRank=frontRank+1;
   PF(frontRank).R=nonPF;
end % end of PF(frontRank)

% Now carry out the sorting of ranks and then update 
[~,frontRanks_Index]=sort(x(:, Rcol));
Sorted_frontRank=x(frontRanks_Index,:); 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Evaluate the crowding distances for each solution for each frontRank   %
% That is, all the non-domonated solutions on the Pareto Front.  %%%%%%%%%%
Qi=0;      % Initialize a counter
for frontRank=1:(length(PF)-1), 
    % Define/initialize a generalized distance matrix 
    dc = [];    past_Q=Qi+1;
    for i=1:length(PF(frontRank).R),
        dc(i,:)=Sorted_frontRank(Qi+i,:);
    end
    Qi=Qi+i;
    % Solutions are sorted according to their fitness/objective values
    fobj_sorted=[];
    for i=1:m, 
        [~, f_Rank]=sort(dc(:,ndim+i));
        fobj_sorted=dc(f_Rank,:);
        % Find the max and min of the fobj values   
        fobj_max=fobj_sorted(length(f_Rank), ndim+i);
        fobj_min=fobj_sorted(1, ndim+i);
        % Calculate the range of the fobj
        f_range=fobj_max-fobj_min;
        % If the solution is at the end/edge, set its distance as infinity
        dc(f_Rank(length(f_Rank)), Rcol+i)=Inf;
        dc(f_Rank(1), Rcol+i) = Inf;
        for j=2:length(f_Rank)-1, 
            fobj2=fobj_sorted(j+1,ndim + i);
            fobj1=fobj_sorted(j-1,ndim + i);  
            % Check the range or special cases
            if (f_range==0),
                dc(f_Rank(j), Rcol+i)=Inf;
            else
            % Scale the range for distance normalization     
            dc(f_Rank(j),Rcol+i)=(fobj2-fobj1)/f_range;
            end
        end % end of j
    end % end of i
    
    % Calculate and update the crowding distances on the Pareto Front
    dist = []; dist(:,1)=zeros(length(PF(frontRank).R),1);
    for i=1:m, 
        dist(:,1)=dist(:,1)+dc(:, Rcol+i);
    end
    % Store the crowding distrance (dc) in the column of Rcol+1=ndim+m+2
    dc(:, Rcol+1)=dist;  dc=dc(:,1:Rcol+1);
    % Update for the output
    xy(past_Q:Qi,:)=dc;  
end  % end of all ranks search/update
sorted_x=xy();    % Output the sorted solutions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% end of non-dominated sorting %%%%%%%%%%%%%%
end

完整代码