MATLAB|模拟植物生长算法（一）

模拟植物生长算法是依据植物向光性原理而设计的智能优化算法，最初是用来解决非线性整数规划问题的，它主要是以植物的生长规则和植物向光性理论为基础建立的概率生长模型，两者结合而形成的优化模式。由于对参数的确定极为简单和宽松，因而具有良好的应用和推广前景，目前主要应用在网络布局优化、设施选址等工程技术领域。（对于参数基本不需要任何设定，只需要慢慢跑就行了，因此它对于初始点的选择有着较为严格的要求，好的初始点会大大缩短其搜索过程）

本文利用一个传统的物流设施选址问题进行引入：

传统的物流设施选址多用重心法以及迭代重心法进行（迭代重心法更为准确，可见论文[1]鲁晓春,詹荷生.关于配送中心重心法选址的研究[J].北方交通大学学报,2000(06):108-110.）我们利用迭代重心法（或者叫精确重心法）得到的结果为（8.39，7.86）。

以此为例对于模拟植物生长算法进行阐述。算法重点内容由以下五部分构成：

（1）植物的生长点（随机初始可行解）

模拟植物生长算法对初始解的要求较低，并且最终的计算结果不依赖于初始可行解，所以可以随机产生一个或者多个初始可行解，并将初始可行解作为植物的初始生长点。（本文用传统的重心法生成初始解，用迭代中心法的结果进行验证）相当于给种子划定一个点，从这个点向四周进行生长。

（2）植物向光性（评价函数）

产生可行解后，就要对可行解进行评价，一般要根据目标函数的特性和植物的向光性特点选择评价函数。植物生长与形态素浓度有关，形态素具有背光的特点，形态素浓度越高，植物生长越快。物流设施选址是以成本最小化为目标函数，恰好与形态素浓度函数相反。因此，论文采用目标函数的负数作为评价函数。

相当于植物接受光照的强弱，接受的光照越强，越有可能向该点生长。

（3）植物生长方向（L系统或T系统）

20世纪60年代末，以简单的重写规则和分枝规则为基础，建立了植物的描述、分析和发育模拟的形式语法，称为L系统。对植物生长作形式化描述，可以根据以下几点进行：1)破土而出的莲杆在一些叫做节的部位长出新枝；大多数新枝上又长出更新的枝，这种分枝行为反复进行；2)不同的枝彼此有相似性，整个植物有自相似结构。系统的简要情形如下：设分枝发生在二维平面上，每次分枝长出的均为单位长，或者旋转一定的角度，如取45°,从种子a开始，采取重写分枝的植物生长规则。T系统与L系统生长情形类似，唯一的不同就是取90°。论文为了简化计算，将采用T系统的方式进行生长。（90度方便计算，但是不够精确，普通模拟植物生长算法固定步长，对于解的搜索不够精确）

（4）形态素浓度（选择和接受解的准则）

植物的向光性问题涉及生物学理论中的形态发生模型，最初的动力学模型是拉什夫斯基、图林等人提出来的叶序模型。叶序问题是关于叶芽细胞、分枝细胞和其他导致叶芽和分枝的分化细胞的生长模式。每个细胞中有均匀的化学物质，其中的一种化学物质是生长激素，称为形态素。这种形态素的浓度随着参量在0和1之间变动，并决定细胞的生长函数是否开始起作用，即细胞分裂，枝芽开始出现。

新的生长点（细胞）产生后，形态素浓度将根据新系统所在环境的改变，重新进行分配。形态素浓度的和是一定的（设为1，且形态素浓度越大，生长的可能性就越大）。为了简化计算，本文将所有现有可生长点的评价函数加总并分配，实际的模拟植物生长算法形态素浓度计算步骤如下：

（5）新枝停止生长（终止准则）

一般设定最大生长次数作为终止准则

算法的流程图如下：

 有了大致的理论，便可以着手MATLAB代码的编写：

clear;
clc;

%% 设置初始参数
t=1;%%计数器
ori=[0,0];%%初始点
bes=ori;
tmax=3000;%%最大生长次数
xi=[0];%%储存已生成的节点x坐标
yi=[0];%%储存已生成的节点y坐标
long=0.3;%%植物茎生长长度

%% 初步计算生长浓度
xi=[xi,1,0,1];
yi=[yi,1,1,0];
xo=xi;
yo=yi;
len=length(xi);
for i =1:len
    g(i)=gm(0,0)-gm(xi(i),yi(i));
end
%%%更新最优值
bes=max(g);
besidx=find(g==bes);
xbes=xi(besidx);

line([xi(1),xi(2)],[yi(1),yi(2)]);
hold on;
line([xi(1),xi(3)],[yi(1),yi(3)]);
hold on;
line([xi(1),xi(4)],[yi(1),yi(4)]);
hold on;

xi(1)=[];
yi(1)=[];%%移除生长点
g(1)=[];
len=length(g);

%%%计算总适应度
tot=0;
for i=1:len
   tot=tot+g(i);
end

%%%计算各点生长素浓度
for i=1:len
   ret(i)=g(i)/tot;
end

%%%随机数选择生长点
rand_num=rand(1);

retn(1)=ret(1);
for i=2:len
   retn(i)=retn(i-1)+ret(i);
end

for i=2:len
    if rand_num<retn(1)
        new=1;
    end
    if rand_num>retn(i-1) && rand_num<retn(i)    
        new=i;
    end
    
end

这一步的目的相当于正式迭代前的预演，展示的是由初始点（0，0）向周围生长的过程，生长了三个点，加入集合“待生长点”，在这三个点中选择下一次生长的点，并将已经生长过的点加入集合“已生长点”，已生长点表明不会有分支向这个点生长。

接着开始正式的迭代操作：

%% 迭代计算
for t=1:tmax
    [xi,yi,xo,yo]=growth(xi,xo,yi,yo,new,long);
    xi(:,new)=[];
    yi(:,new)=[];
    pause(0.0001)
    
    %%%更新总生长素浓度
    len=length(xi);
    for i =1:len
        g(i)=gm(0,0)-gm(xi(i),yi(i));
    end
    
    besg=max(g);
    if besg>bes
        bes=besg;
        besidx=find(g==bes);
        xbes=xi(besidx);
        ybes=yi(besidx);
    end
        
    %%%计算总适应度
    tot=0;
    for i=1:len
        tot=tot+g(i);
    end
    
    %%%计算各点生长素浓度
    for i=1:len
        ret(i)=g(i)/tot;
    end

    %%%随机数选择生长点
    rand_num=rand(1);

    retn(1)=ret(1);
    for i=2:len
       retn(i)=retn(i-1)+ret(i);
    end

    for i=2:len
        if rand_num<retn(1)
            new=1;
        end
        if rand_num>retn(i-1) && rand_num<retn(i)    
            new=i;
        end
    end
    t=t+1;
end
%%输出最终结果
disp(bes);
disp(xbes);
disp(ybes);

 可以发现，这一步与上一步的逻辑十分相似，因此在对于代码熟练的情况下，可以直接从这一步开始编写，减少工作量。

%% 函数，计算适应度,可以根据实际需求进行调整
function [gmi]=gm(a,b)

xz=[4,3,7,13,7,3,11,21];
yz=[14,9,11,8,6,4,4,12];
wz=[10,5,15,5,20,15,10,80];
qz=[5000,5000,5000,5000,5000,5000,5000,2000];
gmi=0;
for i=1:8
    gmi=gmi+sqrt((xz(i)-a).^2+(yz(i)-b).^2).*wz(i).*qz(i);
end

end

这步可以根据需要来调整适应度函数。

%% 函数，设置生长模型，采用L型生长
function [xi,yi,xo,yo]=growth(xi,xo,yi,yo,new,long)
%%%new为新选出的生长点
new1=[xi(new),yi(new)+long];
new2=[xi(new),yi(new)-long];
new3=[xi(new)+long,yi(new)];
new4=[xi(new)-long,yi(new)];
if(sel(xo,yo,new1)==1)
    xi=[xi,new1(1)];
    yi=[yi,new1(2)];
    xo=[xo,new1(1)];
    yo=[yo,new1(2)];
    line([xi(new),new1(1)],[yi(new),new1(2)]);
    hold on;
end
if(sel(xo,yo,new2)==1)
    xi=[xi,new2(1)];
    yi=[yi,new2(2)];
    xo=[xo,new2(1)];
    yo=[yo,new2(2)];
    line([xi(new),new2(1)],[yi(new),new2(2)]);
    hold on;
end
if(sel(xo,yo,new3)==1)
    xi=[xi,new3(1)];
    yi=[yi,new3(2)];
    xo=[xo,new3(1)];
    yo=[yo,new3(2)];
    line([xi(new),new3(1)],[yi(new),new3(2)]);
    hold on;
end
if(sel(xo,yo,new4)==1)
    xi=[xi,new4(1)];
    yi=[yi,new4(2)];
    xo=[xo,new4(1)];
    yo=[yo,new4(2)];
    line([xi(new),new4(1)],[yi(new),new4(2)]);
    hold on;
end

end

这个是L型生长的逻辑，由于本人能力有限，写的...有点儿一言难尽，但是逻辑一看就非常清晰（狗头）

判断某点是否已被生长：

%% 函数，判断新生长点是否已有
function [key]=sel(xi,yi,new)
key=1;%%为1代表为新点，为0代表已有
index1=find(xi==new(1));
index2=find(yi==new(2));
index3=intersect(index1,index2);

if isempty(index3)==0
    key=0;
end
end

最终运行结果如下：

但是模拟植物生长算法仍旧具有许多缺点：

（1）较大的生长空间导致优化效率降低。当生长空间(可行域空间) 较大时，随着生长次数的增大，可生长点集合也随之不断增大，虽然按照植物向光性机理，较优的可生长点仍具有相对较大的概率得到生长的机会，但由于原始的 PGSA 中仅剔除比初始函数值更劣的新增可生长点，较多劣质的生长点会保留在可生长点集合中，大大降低了其优化效率。

（2）初始生长点的选择不当导致收敛于局部最优解。初始生长点不仅影响算法的起点，根据 PGSA的形态素浓度公式，初始生长点的函数值在整个生长过程中都是很重要的参数，影响着算法的计算效率; 且由于PGSA的单一步长搜索机制，需要从初始生长点出发，然后逐步地向着光源渐进，因此当初始生长点在局部最优解附近时，易陷入并收敛于局部最优解。

（3）算法缺乏有效的终止判断机制。原始PGSA是以植物长满整个生长空间(即没有可生长点)为终止条件的，算法往往已经得到了最优解以后仍继续运行，而且实际应用中由于生长空间较大，植物长满整个生长空间所需的计算代价是非常巨大的．

具体改进方法可以参考论文(石开荣,潘文智,姜正荣等.模拟植物生长算法的结构优化新机制[J].华南理工大学学报(自然科学版),2019,47(07):40-48+57.)

本期分享到此结束，由于本人能力有限，有许多缺漏之处，在未来随着我的学习，会逐步增加新的东西，也会持续分享学习的进展。望诸位不吝赐教。