黏菌优化算法SMA在电力系统中虽然有所应用，但是发表的文章还是比较少的，所有这个算法的应用空间还很广，推荐推荐：
对于这两篇文献，我写过相关的文章，先回顾一下，然后再开始今天的话题：黏菌优化算法SMA 。
求解热电联产经济调度问题的改进遗传与粒子群算法 |
Simulink|电动汽车、永磁电动机建模与仿真 |
多目标粘液霉菌算法（Matlab代码实现） |
目录
1 黏菌模型——人类成功模拟整个宇宙网络！
2 黏液霉菌优化算法数学模型
2.1 接近食物
2.2 包裹食物
2.3 抢食物
3 算法流程
4 Matlab代码实现
4.1 Matlab代码
4.2 运行结果
1 黏菌模型——人类成功模拟整个宇宙网络！
本部分参考：百家号 #科学了不起# 系列征文赛：计算机模拟单细胞生物模型,科学家沸腾了！利用生物觅食模型，人类成功模拟整个宇宙网络！...
好消息，地球上的黏菌模型成功帮助天文学家绘制了连接整个宇宙星系的宇宙网络！
黏菌，也称多头绒泡菌，是一种单细胞生物，它通过建立复杂的丝状网络来寻找食物。加州大学圣克鲁斯分校的研究人员受黏液霉菌生长模式的启发，利用计算机模型，追踪了星系之间延伸数光年的网状细丝。
来自加州大学圣克鲁兹分校的这项研究的主要作者乔·伯切特(Joe Burchett)：“黏液霉菌创造了一个优化的运输网络，找到了连接食物来源的最有效途径，同样的，在宇宙网络中，结构的增长产生的网络在某种意义上也是最优的。其基本过程是不同的，但它们产生了类似的数学结构。”
为了创建新的模型，该团队使用了来自斯隆数字巡天计划(Sloan Digital Sky Survey)的数据和柏林艺术家塞奇·詹森(Sage Jenson)的作品，他的艺术可视化基于一种模拟黏液霉菌生长的算法。根据声明，研究人员将这种新算法命名为蒙特卡罗绒泡机。
宇宙中的物质分布在由巨大空洞分隔的星系间细丝构成的网状网络中，而星系形成于这些细丝相交的地方，恰恰是物质最集中的地方。然而，这些在星系之间延伸的细丝基本上是看不见得，因为它们是由暗物质组成的，这是一种不发光或不释放能量的物质，但却占了宇宙质量的大约85%。
研究人员用来自博尔修伊-普朗克宇宙学模拟的数据对新算法进行了测试。这个模拟是由加州大学圣克鲁兹分校的物理学教授JoelPrimack开发的，它被用来模拟暗物质的“晕”，即星系形成的地方以及连接整个宇宙星系的细丝。结果表明，新的黏液霉菌算法的结果与暗物质模拟的结果紧密一致。
黏液霉菌绒泡菌在探索食物环境时，形成了一个相互连接的管道网络。科学家受其生长模式启发的算法看到了连接所有星系的宇宙网络的结构。
“从45万个暗物质晕开始，我们可以在宇宙学模拟中得到一个近乎完美的密度场，”该研究的合著者、加州大学圣克鲁兹分校(UC Santa Cruz)计算媒体博士后研究员奥斯卡·埃莱克(Oskar Elek)在声明中说。
研究人员还使用了来自哈勃太空望远镜的宇宙起源光谱仪的数据，该光谱仪用于研究吸收或发射光的物体。根据这份声明，星系间的气体在穿过它的光谱中留下了独特的吸收特征。因此，这些哈勃数据揭示了星系间空间的气体特征。根据声明，这些气体信号在靠近细丝中心的地方更强烈，在那里密集的物质堆积形成了新的星系。
随后研究人员利用斯隆数字巡天计划(SloanDigitalSkySurvey)和哈勃太空望远镜(HubbleSpaceTelescope)的数据，使用一种算法来模拟地球上黏液霉菌的生长模式，以便更好地模拟宇宙网络的结构。
伯切特在声明中说：“现在，我们第一次可以量化星系间介质的密度，从遥远的宇宙网络细丝的边缘到炽热、稠密的星系团内部。”“这些结果不仅证实了宇宙学模型所预测的宇宙网络的结构，它们还为我们提供了一种方法，通过将其与构成星系的气藏联系起来，来提高我们对星系演化的理解。”
因此，基于黏液模型的新算法帮助了天文学家在更大的尺度上建立可视化宇宙网络。
2 黏液霉菌优化算法数学模型
2.1 接近食物
将粘菌的接近行为建模为数学方程，提出以下规则来模拟收缩模式：
其中
是一个范围为
的参数，
从 1 线性减小到 0。
代表当前迭代，
代表当前发现的气味浓度最高的个体位置，
代表粘菌的位置，
和
代表从群体中随机选择的两个个体，
代表粘菌的重量。
的公式如下：
其中
,表示
的适应度，
表示在所有迭代中获得的最佳适应度。
的公式如下:
的公式如下所示:
2.2 包裹食物
更新黏菌位置的数学公式如下:
其中
和
表示搜索范围的上下边界，表示[0,1]中的随机值。
2.3 抢食物
的值在
之间随机振荡，并随着迭代次数的增加逐渐接近零。
值在 [-1,1] 之间振荡，最终趋于零.
Matlab代码表述为：
%=============欢迎关注公众号：电力系统与算法之美=======
% max _ iter：最大迭代次数，N：种群大小，收敛曲线：收敛曲线
function [Destination_fitness,bestPositions,Convergence_curve]=SMA(N,Max_iter,lb,ub,dim,fobj)
%% 初始化位置
bestPositions=zeros(1,dim);
Destination_fitness=inf;%将此更改为 -inf 以解决最大化问题
AllFitness = inf*ones(N,1);%记录所有粘菌的适应度
weight = ones(N,dim);%每个粘菌的适应度权重
%% 初始化随机解集
X=initialization(N,dim,ub,lb);
Convergence_curve=zeros(1,Max_iter);
it=1; %迭代次数
lb=ones(1,dim).*lb; % 变量下限
ub=ones(1,dim).*ub; % 变量上限
z=0.03; % 参数
%% 主循环
while it <= Max_iter
%=====适应度排序======
for i=1:N
% 检查解决方案是否超出搜索空间并将其带回
Flag4ub=X(i,:)>ub;
Flag4lb=X(i,:)<lb;
X(i,:)=(X(i,:).*(~(Flag4ub+Flag4lb)))+ub.*Flag4ub+lb.*Flag4lb;
AllFitness(i) = fobj(X(i,:));
end
[SmellOrder,SmellIndex] = sort(AllFitness);
worstFitness = SmellOrder(N);
bestFitness = SmellOrder(1);
S=bestFitness-worstFitness+eps; %加上 eps 以避免分母为零
%====计算每个粘菌的适应度权重=====
for i=1:N
for j=1:dim
if i<=(N/2)
weight(SmellIndex(i),j) = 1+rand()*log10((bestFitness-SmellOrder(i))/(S)+1);
else
weight(SmellIndex(i),j) = 1-rand()*log10((bestFitness-SmellOrder(i))/(S)+1);
end
end
end
%====更新最佳适应度值和最佳位置=====
if bestFitness < Destination_fitness
bestPositions=X(SmellIndex(1),:);
Destination_fitness = bestFitness;
end
a = atanh(-(it/Max_iter)+1);
b = 1-it/Max_iter;
%====更新搜索代理的位置=====
for i=1:N
if rand<z
X(i,:) = (ub-lb)*rand+lb;
else
p =tanh(abs(AllFitness(i)-Destination_fitness));
vb = unifrnd(-a,a,1,dim);
vc = unifrnd(-b,b,1,dim);
for j=1:dim
r = rand();
A = randi([1,N]); % 从总体中随机选择两个位置
B = randi([1,N]);
if r<p
X(i,j) = bestPositions(j)+ vb(j)*(weight(i,j)*X(A,j)-X(B,j));
else
X(i,j) = vc(j)*X(i,j);
end
end
end
end
Convergence_curve(it)=Destination_fitness;
it=it+1;
end
end
3 算法流程
算法流程：
Step1.初始化种群，设定相应参数。
Step2.计算适应度值，并且排序。
Step3.利用式（1），更新种群位置。
Step4.计算适应度值，并且更新全家最优位置，当前最优位置。
Step5.是否达到结束条件，如果达到则输出最优结果，否则重复执行步骤2-5.
4 Matlab代码实现
4.1 Matlab代码
%=============欢迎关注公众号：电力系统与算法之美=======
clear all
close all
clc
N=30;
Function_name='F2'; % 目标函数的名称，范围从 F1-F13（也可以换成自己的目标函数）
T=500; %最大迭代次数
dimSize = 30; %维度
%=======调用目标函数===============
[lb,ub,dim,fobj]=Get_Functions_SMA(Function_name,dimSize);
[Destination_fitness,bestPositions,Convergence_curve]=SMA(N,T,lb,ub,dim,fobj);
%========绘制目标空间======
figure,
hold on
semilogy(Convergence_curve,'Color','b','LineWidth',4);
title('收敛曲线')
xlabel('迭代次数');
ylabel('最优解');
axis tight
grid off
box on
legend('SMA')
display(['SMA的最佳位置是: ', num2str(bestPositions)]);
display(['SMA的最佳适应度是: ', num2str(Destination_fitness)]);
%=============欢迎关注公众号：电力系统与算法之美=======
% max _ iter：最大迭代次数，N：种群大小，收敛曲线：收敛曲线
function [Destination_fitness,bestPositions,Convergence_curve]=SMA(N,Max_iter,lb,ub,dim,fobj)
%% 初始化位置
bestPositions=zeros(1,dim);
Destination_fitness=inf;%将此更改为 -inf 以解决最大化问题
AllFitness = inf*ones(N,1);%记录所有粘菌的适应度
weight = ones(N,dim);%每个粘菌的适应度权重
%% 初始化随机解集
X=initialization(N,dim,ub,lb);
Convergence_curve=zeros(1,Max_iter);
it=1; %迭代次数
lb=ones(1,dim).*lb; % 变量下限
ub=ones(1,dim).*ub; % 变量上限
z=0.03; % 参数
%% 主循环
while it <= Max_iter
%=====适应度排序======
for i=1:N
% 检查解决方案是否超出搜索空间并将其带回
Flag4ub=X(i,:)>ub;
Flag4lb=X(i,:)<lb;
X(i,:)=(X(i,:).*(~(Flag4ub+Flag4lb)))+ub.*Flag4ub+lb.*Flag4lb;
AllFitness(i) = fobj(X(i,:));
end
[SmellOrder,SmellIndex] = sort(AllFitness);
worstFitness = SmellOrder(N);
bestFitness = SmellOrder(1);
S=bestFitness-worstFitness+eps; %加上 eps 以避免分母为零
%====计算每个粘菌的适应度权重=====
for i=1:N
for j=1:dim
if i<=(N/2)
weight(SmellIndex(i),j) = 1+rand()*log10((bestFitness-SmellOrder(i))/(S)+1);
else
weight(SmellIndex(i),j) = 1-rand()*log10((bestFitness-SmellOrder(i))/(S)+1);
end
end
end
%====更新最佳适应度值和最佳位置=====
if bestFitness < Destination_fitness
bestPositions=X(SmellIndex(1),:);
Destination_fitness = bestFitness;
end
a = atanh(-(it/Max_iter)+1);
b = 1-it/Max_iter;
%====更新搜索代理的位置=====
for i=1:N
if rand<z
X(i,:) = (ub-lb)*rand+lb;
else
p =tanh(abs(AllFitness(i)-Destination_fitness));
vb = unifrnd(-a,a,1,dim);
vc = unifrnd(-b,b,1,dim);
for j=1:dim
r = rand();
A = randi([1,N]); % 从总体中随机选择两个位置
B = randi([1,N]);
if r<p
X(i,j) = bestPositions(j)+ vb(j)*(weight(i,j)*X(A,j)-X(B,j));
else
X(i,j) = vc(j)*X(i,j);
end
end
end
end
Convergence_curve(it)=Destination_fitness;
it=it+1;
end
end
%=============欢迎关注公众号：电力系统与算法之美=======
% 初始化
function Positions=initialization(SearchAgents_no,dim,ub,lb)
Boundary_no= size(ub,2); % 边界数
%如果所有变量的边界相等，对ub和lb都输入一个标志号
if Boundary_no==1
Positions=rand(SearchAgents_no,dim).*(ub-lb)+lb;
end
%如果每个变量有不同的 lb 和 ub
if Boundary_no>1
for i=1:dim
ub_i=ub(i);
lb_i=lb(i);
Positions(:,i)=rand(SearchAgents_no,1).*(ub_i-lb_i)+lb_i;
end
end
%=============欢迎关注公众号：电力系统与算法之美=======
%=====目标函数==============
function [lb,ub,dim,fobj] = Get_Functions_SMA(F,DimValue)
switch F
case 'F1'
fobj = @F1;
lb=-100;
ub=100;
dim=DimValue;
case 'F2'
fobj = @F2;
lb=-10;
ub=10;
dim=DimValue;
case 'F3'
fobj = @F3;
lb=-100;
ub=100;
dim=DimValue;
case 'F4'
fobj = @F4;
lb=-100;
ub=100;
dim=DimValue;
case 'F5'
fobj = @F5;
lb=-30;
ub=30;
dim=DimValue;
case 'F6'
fobj = @F6;
lb=-100;
ub=100;
dim=DimValue;
case 'F7'
fobj = @F7;
lb=-1.28;
ub=1.28;
dim=DimValue;
case 'F8'
fobj = @F8;
lb=-500;
ub=500;
dim=DimValue;
case 'F9'
fobj = @F9;
lb=-5.12;
ub=5.12;
dim=DimValue;
case 'F10'
fobj = @F10;
lb=-32;
ub=32;
dim=DimValue;
case 'F11'
fobj = @F11;
lb=-600;
ub=600;
dim=DimValue;
case 'F12'
fobj = @F12;
lb=-50;
ub=50;
dim=DimValue;
case 'F13'
fobj = @F13;
lb=-50;
ub=50;
dim=DimValue;
end
end
% F1
function o = F1(x)
o=sum(x.^2);
end
% F2
function o = F2(x)
o=sum(abs(x))+prod(abs(x));
end
% F3
function o = F3(x)
dim=size(x,2);
o=0;
for i=1:dim
o=o+sum(x(1:i))^2;
end
end
% F4
function o = F4(x)
o=max(abs(x));
end
% F5
function o = F5(x)
dim=size(x,2);
o=sum(100*(x(2:dim)-(x(1:dim-1).^2)).^2+(x(1:dim-1)-1).^2);
end
% F6
function o = F6(x)
o=sum(abs((x+.5)).^2);
end
% F7
function o = F7(x)
dim=size(x,2);
o=sum([1:dim].*(x.^4))+rand;
end
% F8
function o = F8(x)
o=sum(-x.*sin(sqrt(abs(x))));
end
% F9
function o = F9(x)
dim=size(x,2);
o=sum(x.^2-10*cos(2*pi.*x))+10*dim;
end
% F10
function o = F10(x)
dim=size(x,2);
o=-20*exp(-.2*sqrt(sum(x.^2)/dim))-exp(sum(cos(2*pi.*x))/dim)+20+exp(1);
end
% F11
function o = F11(x)
dim=size(x,2);
o=sum(x.^2)/4000-prod(cos(x./sqrt([1:dim])))+1;
end
% F12
function o = F12(x)
dim=size(x,2);
o=(pi/dim)*(10*((sin(pi*(1+(x(1)+1)/4)))^2)+sum((((x(1:dim-1)+1)./4).^2).*...
(1+10.*((sin(pi.*(1+(x(2:dim)+1)./4)))).^2))+((x(dim)+1)/4)^2)+sum(Ufun(x,10,100,4));
end
% F13
function o = F13(x)
dim=size(x,2);
o=.1*((sin(3*pi*x(1)))^2+sum((x(1:dim-1)-1).^2.*(1+(sin(3.*pi.*x(2:dim))).^2))+...
((x(dim)-1)^2)*(1+(sin(2*pi*x(dim)))^2))+sum(Ufun(x,5,100,4));
end
function o=Ufun(x,a,k,m)
o=k.*((x-a).^m).*(x>a)+k.*((-x-a).^m).*(x<(-a));
end
4.2 运行结果
